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ABSTRACT 


Composite  reliability  is  strongly  dependent  on  the  fiber  strength 
distribution.  Current  methods  of  gathering  statistics  through  single  fiber  failure 
methods  are  inefficient  and  costly.  This  thesis  develops  a  testing  method  from 
which  the  fiber  statistics  in  the  form  of  Weibull  parameters  can  be  accurately 
extracted  from  bundle  failure  tests. 

The  values  obtained  from  the  bundle  experiment  as  compared  to  known 
single  fiber  test  bench  mark  parameters  were  practically  indistinguishable .  The 
confidence  of  the  results  stems  from  a  thorough  analysis  of  the  bundle 
mechanisms  and  minimized  contaminates  which  can  disturb  the  strength 
distribution. 
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I.  INTRODUCTION 


A  composite,  by  definition,  is  any  two  or  more  materials  combined  on  a 
macroscopic  scale  to  form  a  useful  material  [Ref  1J.  Today,  in  one  form  or 
another,  composites  are  being  used  on  every  level  of  our  society.  This  thesis 
will  focus  on  a  class  of  composites  known  as  "fibrous  composites,"  today  these 
are  typically  used  in  high  tech  applications,  often  high  performance  military 
aircraft.  Fiber  composites  consist  of  very  thin  fibers  with  high  aspect  ratios  and 
high  strengths  imbedded  in  a  matrix,  commonly  a  polymeric  matrix  such  as 
epoxy.  Glass/epoxy  is  an  example  of  a  frequently  utilized  composite  mixture 
that  has  been  in  use  for  many  years. 

One  of  the  many  advantages  this  class  of  material  provides  is  high 
strength-to-weight,  though  due  to  a  lack  of  thorough  understanding  in  failure 
mechanics,  composites  have  been  limited  to  secondary  load  carrying  structural 
designs  such  as  skins  and  control  surfaces.  In  recent  years  (last  20),  composites 
are  expanding  into  a  much  wider  role. 

A.  REQUIREMENTS  IN  ADVANCED  STRUCTURES 

Since  the  introduction  of  glass/epoxy,  advances  in  production  techniques  and 
developments  in  the  textile  industry  have  resulted  in  fiber  materials  with  much 
improved  strengths  and  stiffness.  Two  of  the  early  break-throughs  were 
graphite  and  boron  fibers  which  have  typical  strengths  and  stiffness  much 
higher  than  convention  ductile  isotropic  metals.  These  fibers  allow  for  new 


design  limits  but  also  give  rise  to  the  requirement  of  new  design  technique 
utilizing  a  material  of  an  unconventional  material  redundancy  with  anisotropy. 

Employing  the  fiber/matrix  as  unidirectional  lamina  to  form  multidirectional 
laminate  gave  the  designer  the  tools  to  make  structural  designs  subjected  to 
multiple  loading  conditions  a  manageable  task.  This  idea  shifted  the  attention 
of  the  designer  from  the  microscopic  level  to  the  macroscopic  level. 

Utilizing  the  high  strength  and  stiffness,  fibers  application  for  structural 
enhancement  were  incorporated  in  many  of  the  military  designs  in  the  early 
70’s  such  as  the  F-lll  afl  fuselage  structure.  Testing  of  these  designs  were 
accomplished  by  component.  This  method  is  time  consuming  and  expensive  but 
must  be  done.  To  gain  any  type  of  reliability  assurance  many  tests  are 
required.  It  soon  becomes  apparent  that  testing  costs  could  out  pace  production 
costs. 

B.  FUTURE  DESIGN  REQUIREMENTS 

As  weaponry  becomes  more  and  more  sophisticated  the  demands  on 
composite  materials  for  a  wide  variety  of  uses  are  ever  increasing.  New  designs 
often  are  only  limited  by  the  structural  limitations  of  the  materials  in  use. 
Such  restrictions  might  include  the  thickness  of  a  wing  or  stabilizer,  the  sweep 
of  the  wing,  or  the  "g"  loading.  Lately  new  design  specifications  have  included 
the  use  of  composites  for  the  radiation  energy  absorbing  properties.  Though  the 
Navy’s  interest  in  composites  spans  a  multitude  of  uses,  one  special  interest  is 
in  the  area  of  structural  enhancement,  such  as  rocket  motor  cases,  pressure 
vessels  for  submarine  flasks  and  jet  aircraft  pilot  ejection  seat. 
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In  the  last  few  decades  composites  have  advanced  dramatically,  and  are 
beginning  to  provide  improved  alternatives  to  more  conventional  methods  of 
design  normally  utilizing  isotropic  alloys.  Fiber  reinforced  plastics  have 
attracted  a  large  amount  of  attention.  Improved  techniques  in  the  textile 
industry  have  lead  to  more  uniform  fibers  from  bases  of  precursors  for  high 
performance  graphite’s  with  strengths  and  stiffness  exceeding  high  yield  steels. 

As  utilization  goes  up  the  increased  use  of  new  fibers  in  composites  and  the 
cost  of  system  failure  testing  demand  better  methods  for  reliability  evaluation. 
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II.  BACKGROUND 


A.  COMPOSITE  LOAD  SHARING 

A  brittle  material  generally  has  a  higher  strength  to  density  ratio  but  a 
lack  of  ductility  limits  the  structural  uses.  Ductility  resists  crack  propagation 
from  a  region  of  damage.  Metals  have  a  high  ductility  and  do  provide  practical 
strength  and  stiffness  properties  for  many  structural  applications.  Fibers  do  not 
use  ductility  to  resist  failure  [Ref.  2:pp.  1-19],  because  ductility  is  not  required 
the  fiber  can  be  synthesized  from  brittle  materials. 

The  maximum  theoretical  elastic  strain  for  carbon  hydrogen  bonds  is 
approximately  £*„  =  0.1,  which  is  one  order  of  magnitude  greater  than  1 

present  day  strong  elastic  reinforcing  fibers. 

Presently  fibers  can  not  obtain  strengths  any  where  near  the  theoretical 
ideal  strength  of  the  material  but  are  an  order  of  magnitude  greater  than  ! 

conventional  bulk  materials.  Fiber  strength  can  also  be  attributed  to  the  small 
fiber  diameter,  which  limits  the  flaw  sizes. 

Composites  are  geometrically  advantageous  to  homogeneous  material  such  as  ! 

an  all  metal  or  all  ceramic  material,  for  its  ability  to  utilize  the  high  strengths 
and  stiffness  of  brittle  materials  while  maintaining  strength  redundancy  through 
the  matrix.  ( 

The  strength  interaction  of  the  fiber  and  matrix  is  the  key  to  reliability 
effectiveness  of  a  composite.  Consider  several  fibers  in  parallel  without  matrix. 

The  load  is  equally  shared  by  each  fiber.  If  one  fiber  should  fail  the  failed  fiber  j 

is  completely  ineffective  and  no  longer  supports  any  load.  The  load  of  the 


I 


i 

remaining  fibers  increases  equally.  With  a  matrix  present  (Figure  1)  a  crucial 
phenomenon  occurs.  At  the  point  of  fracture  the  more  ductile  matrix  absorbs 
the  load  through  shear  stresses  along  adjacent  fibers  and  back  to  the  broken  I 

fiber  itself.  Near  the  fracture  site  high  shear  stress  develope  in  the  matrix 
along  the  broken  fiber  but  quickly 
dissipate  as  the  normal  stress 
increases.  On  the  two  surrounding 
fibers  a  noticeable  stress  increase  is 
observed  adjacent  to  the  fracture. 

The  short  distance  from  the  broken 
end  of  the  fiber  to  the  point  that 
fiber  is  carrying  a  full  stress  is  called 
the  ineffective  length  [Ref.  3] 

The  ability  of  the  composite  to 
adjust  to  multiple  fiber  failures  is  it’s 
strength  redundancy,  and  is 
dependent  on  the  characteristics  of 

both  the  matrix  and  fiber.  There  are  several  fracture  mechanisms  all  of  which 
are  a  function  of  the  bond  strength  (shear  strength  of  interface),  the  matrix 
shear  strength  and  the  distribution  of  the  fiber  strength.  The  most  critical  is 
the  fiber  strength  distribution.  The  distribution  of  the  fiber  determines  the 
composite  reliability. 
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B.  STATISTICS  AND  RELIABILITY 

The  complexities  of  composites  make  the  task  of  analysis  virtually 
impossible.  New  fibers  are  being  developed  every  year  and  attempting  to 
generate  thorough  material  testing  studies  as  has  been  done  with  conventional 
materials  is  not  realistic.  But  composite  reliability  can  be  estimated  through 
probability  studies.  Probability  of  the  composite  must  be  based  or  inferred  from 
accurate  statistics  of  the  fiber  distribution. 

To  determine  the  reliability  of  the  composite,  the  strength  distribution  of 
the  fiber  must  be  known.  The  distribution  can  be  quantified  by  strength 
statistics  (a  Empirical  Probability  Distribution  Function  (EPDF),  or  Empirical 
Cumulative  Distribution  Function  (ECDF)).  For  computational  efficiency  the 
empirical  data  can  then  be  modeled  with  an  analytical  function  and  thereby 
characterize  the  distribution  by  parameters  of  the  function.  To  define  the 
parameters  the  distribution  function  must  be  identified  by  model  type  (i.e., 
Normal,  Weibull,  etc.).  Using  the  model  Probability  plot;  the  closeness  of  fit  of 
the  data  can  be  judge  and  therefore  the  appropriateness  of  the  model.  Selecting 
the  model  allows  for  the  calculation  of  the  model  parameters. 

The  relationship  between  the  fiber  statistics  and  the  composite  probability  is 
represented  in  a  linearized  Probability  plot  (Figure  2). 

Reliability  is  mostly  concerned  with  the  lower  tail  of  the  distribution,  this 
also  being  the  most  difficult  to  define.  Because  of  the  redundancy  provided  by 
the  matrix  induced  load  sharing  the  composite  is  considerably  more  reliable 
than  it’s  constituent  fibers  in  the  lower  tail.  As  the  fiber  statistics  fluctuate  in 
the  extreme  lower  tail  the  composite  reliability  fluctuates  with  it.  This  can  be 
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Probability  Plot 
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Figure  2.  Comparison  of  Fiber  and  Composite  in  a 
Probability  Plot. 

shown  in  a  representation  of  a  probability  plot  in  Figure  3.  There  is  a  critical 
point  where  the  fiber  plot  stabilizes,  the  linear  region  in  this  plot  determines 
the  accurate  statistics  of  the  lower  bound  reliability  (worse  case)  of  the 
composite. 

Fibers  have  a  characteristic  which  helps  to  identify  the  analytical 
distribution  model  that  will  closely  define  the  empirical  statistic.  Traditionally, 
the  weakest  link  model  (Ref.  4]  depicted  the  fiber  as  being  made  of  many  small 
segments  linked  together  much  like  a  chain.  The  segments  each  have  an 
intrinsic  strength,  none  being  exactly  equal.  When  the  fiber  is  stressed  the 
weakest  will  fail  and  the  total  fiber  fails.  The  Weibull  distribution  model 
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defines  this  weakest  link  behavior  and  is  commonly  used  in  modeling  fiber 
strength. 


1.  Weibull  Density  Distribution  Function 

The  Two  parameter  Weibull  cumulative  distribution  function  (CDF)  is 
defined  by  the  following  equation. 


F(x)  =  1  -  exp  (-(x/B)°),  x,  0,  a  >0 


i 


(2.1) 


where:  x  =  random  variable  (stress,  load) 

B  =  scale 

a  =  shape  I 


The  parameter  alpha  defines  the  shape  (spread  and  skew)  of  the 

distribution.  Alpha  less  than  3.5  the 
plot  is  skewed  positively,  greater  than 
3.5  it  is  skewed  negatively  and 
neutral  at  3.5.  The  alpha  parameter 
determines  the  strength  scatter  of  the 
fiber,  an  increasing  alpha  defines  less 
scatter  and  higher  strength 
uniformity.  The  beta  parameter 
defines  the  central  tendency  of  the 
distribution  and  is  relatable  to  the 
mean.  Alpha  for  a  fiber  type  is 
constant  for  the  weak  link  physics, 
whereas  beta  changes  with  the  fiber  length.  An  increase  in  length  causes  a 
decrease  in  beta  (a  weaker  fiber). 
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The  composite  reliability  can  be  inferred  from  the  distribution  of  a  fiber 
strength.  To  gain  high  reliability  with  high  confidence  levels  a  proportionate 
amount  of  statistics  must  be  obtained  (i.e.  for  a  reliability  level  of  10-3  requires 
104  statistics). 

Presently  fiber  statistics  are  gathered  by  single  filament  testing  (Figure 
4).  First  a  single  fiber  specimen  is  prepared,  since  a  fiber  has  such  a  small 
diameter,  (AS-4  diameter  -  10  micron)  it  is  very  fragile.  The  fiber  is  placed  in 
a  cardboard  frame.  Once  the  specimen  is  loaded  in  the  test  apparatus  the 
cardboard  frame  is  severed  by  a 
heated  wire  followed  by  load 
application  and  interpretation  of  the 
load  deformation  data.  This  is  a 
time  consuming  effort  which 
produces  only  one  failure  strength 
statistic.  To  gain  the  number  of 
statistics  necessary  would  take 

thousands  of  tests,  as  an  example  man-safe  reliability  is  10  s  therefore  to 
achieve  the  necessary  statistics  would  require  107  tests. 

A  new  approach  being  explored  is  to  test  many  fibers  in  parallel  (a 
bundle)  and  recover  the  fiber  characteristics  from  the  failure  test.  This  would 
provide  a  large  statistical  base  required  and  could  be  accumulated  with  a  much 
higher  efficiency.  With  a  validated  method  the  analysis  of  fiber  reliability  could 
keep  pace  with  the  development  of  new  fibers  with  reduced  effort. 
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2.  Size  Effect. 

The  extreme  lower  tail  is  difficult  to  determine  accurately  for  the  fiber. 
For  a  Weibull  failure  physics  and  the  resulting  strength  distribution  the  lower 
tail  can  be  estimated  by  assuming  weakest  link  [Ref.  4].  The  parameter  beta 
is  a  function  of  the  fiber  length  and  is  defined  by  equation  2.2 

ii2  =  e.d./lj)*^  (2.2) 

where  B,  =  beta  for  1, 
fi2  =  beta  for  12 
a  s=  a,  = 

A  detailed  development  can  be  seen  in  Appendix  A.  Using  this 
relationship  parameters  for  beta  can  be  retrieved  analytically  where  practical 
testing  is  impossible  and  the  extreme  lower  tail  is  estimated. 

This  thesis  will  attempt  to  measure  the  lower  tail  via  bundle  testing,  as 
opposed  to  the  aforementioned  estimation,  (fiber  statistics  via  single  filament 
testing).  This  investigation  explores  the  methodology  needed  in  testing  a  bundle 
of  fibers  to  failure.  By  recording  the  load  and  displacement  it  is  expected  that 
the  underlying  distribution  of  the  fibers  in  that  bundle  may  be  extracted.  If, 
for  example,  a  bundle  has  3k  of  fibers,  the  distribution  is  based  on  3k  of 
statistics.  As  the  statistical  base  approaches  the  complement  of  the  reliability, 
the  accuracy  of  the  estimation  increases. 


III.  TECHNICAL  ISSUES  AND  RESOLUTIONS 


A  bundle  of  fibers  tested  to  failure  presents  several  experimental  obstacles 
which  effects  the  analysis.  There  are  three  experimental  consequences  to 
testing  a  bundle  in  parallel  that  are  not  experienced  in  testing  a  single  fiber  in 
series.  These  phenomenon  which  may  disturb  the  failure  distribution,  are  1)  a 
protective/adhesive  enhancement  coating  called  sizing  which  is  placed  on  the 
filament  spools  during  production,  2)  friction  which  occurs  during  testing  after  a 
portion  of  the  bundle  has  failed  and  3)  slack  which  results  from  specimen 
fabrication.  Each  of  these  is  discussed  in  the  following  paragraphs. 

A.  SIZING 

Fibers  are  packaged  on  spools  of  bundles,  in  order  to  prevent  the  fibers 
from  entangling  a  small  amount  of  a  chemical  liquid  called  sizing  is  coated  on 
the  bundle.  The  effect  of  this  sizing  on  bundle  testing  is  unknown  and  must  be 
address  at  sometime.  The  testing  in  this  experiment  was  done  leaving  sizing  on 
the  bundle.  A  method  to  remove  the  sizing  chemical  without  damage  to  the 
fibers  needs  to  be  studied  but  is  not  addressed  here. 

B.  FRICTION 

During  the  bundle  test  prior  to  fiber  failure  there  is  negligible  friction,  this  I 

is  due  to  the  equal  strain  between  fibers.  As  the  fibers  fail  the  strain  of  the 
failed  filaments  return  to  zero  leaving  the  loose  ends  entangled  in  the 
remaining  fibers,  this  places  added  strain  on  those  fibers.  The  resulting  effect  I 

shows  up  as  reduced  failure  strengths  and  clustering  of  the  these  failures.  In 
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the  relatively  short  gauge  lengths  friction  is  expected  to  be  small  and  very 
pronounced  in  the  long  gauge  lengths.  The  effects  of  lubricating  the  bundle 
with  oil  during  testing  will  be  observed  to  see  if  the  effect  is  reduced. 

C.  SLACK 

During  specimen  preparation  and  placement  in  the  test  apparatus  the 
alignment  of  the  fibers  becomes  misaligned  with  unequal  length.  The  effective 
result  is  an  initial  nonlinear  load/displacement  which  also  effect  the  failure 
distribution.  Another  concern  is  the  deformation  in  the  actual  load  train  of  the 
test  apparatus.  This  is  referred  to  as  system  compliance.  Compliance  was 
found  to  be  a  function  of  the  load  and  can  be  removed  from  the  resulting  data. 
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IV.  EXPERIMENTATION 


The  experimental  setup  consisted  of  a  Material  Tester  (INSTRON  4206) 
and  an  integrated  data  acquisition  system  consisting  of  an  Instron  Control 
Console,  IBM  PC/AT  and  software,  a  detail  description  and  a  list  of  procedures 
is  in  Appendix  B.  The  object  of  the  experiment  was  to  test  a  current  composite 
fiber  type  (used  in  many  Navy  aircraft.)  with  a  range  of  several  gauge  lengths 
and  attempt  to  extract  the  underlying  distribution  of  the  fiber.  The  sample 
preparation  is  described  in  Appendix  B.  Each  bundle  sample  was  tested,  the 
data  stored,  and  quick  analysis  conducted  to  ensure  no  improper  testing  method 
were  corrupting  the  data. 

A.  BUNDLE  TESTING 

The  fiber  material  tested  was  made  from  a  Hercules  Magnamite  high 
strength  graphite,  type  AS-4  spool  145  in  3k  bundles  (i.e.  there  are  3,000 
filaments  in  a  bundle).  Several  different  gauge  lengths  were  tested  in  order  to 
analyze  the  effects  of  friction,  slack,  compliance  on  the  fiber  size.  Size  affect 
(Appendix  A)  should  be  able  to  correlate  the  location  parameter  beta  between 
fiber  lengths  if  weakest  link  physics  is  applicable.  In  all,  nineteen  bundle 
samples  where  tested.  The  test  procedure  was  conducted  as  specified  in 
Appendix  B.  The  complete  test  log  is  listed  in  Table  I. 

The  bundle  samples  were  handled  with  extreme  care  so  as  not  to  damage 
the  any  of  the  fibers.  Moving  the  bundles  from  assembly  bench  to  the  storage 
bench  and  then  to  the  testing  machine  presented  the  greatest  hazard.  The 
difficulty  came  when  picking  up  the  sample.  This  was  done  manually  and 


Table  I.  TEST  LOG. 


Sample  # 

G.L. 

(cm) 

Oil/ 

Dry 

Comments 

090901 

5.0 

dry 

Good 

090902 

0.5 

dry 

No  failure-slip 

090903 

0.5 

dry 

No  failure-slip 

090904 

0.5 

dry 

No  failure-slip 

090905 

0.5 

dry 

No  failure-slip 

090906 

0.5 

dry 

No  failure-slip 

090907 

5.0 

dry 

Slip 

090908 

5.0 

dry 

Good 

090909 

5.0 

oil 

Good 

090910 

2.5 

dry 

Possible  damage 

100901 

2.5 

dry 

Good 

100902 

2.5 

oil 

Slip 

100903 

2.5 

oil 

Slip 

100904 

50.0 

dry 

Good  -  (friction) 

100905 

50.0 

oil 

Good  -  (friction 
reduced) 

100906 

50.0 

oil 

Possible  damage 

100907 

25.0 

oil 

Good 

100908 

25.0 

dry 

Good 

bending  the  bundles  slightly  was  unavoidable.  The  long  bundles  (500  and  250 
mm)  were  less  vulnerable  to  damage  because  of  the  relatively  large  aspect  ratio 
but  these  were  difficult  to  steady  during  transportation  and  tended  to  bounce 
around. 

Placing  the  sample  in  the  grips  also  presented  operational  difficulties  when 
tightening  the  upper  grip  with  a  socket  wrench.  An  attempt  was  made  to 
manually  counter  torque  the  grip,  but  this  tended  to  be  jerky  at  best.  Also  the 
grip  pressure  was  difficult  to  control  between  samples.  The  alignment  of  the 
bundle  in  the  grips  was  another  concern.  This  was  done  by  sight.  If  not  done 
properly  the  bundle  failure  mechanism  might  be  influenced.  Manually  reducing 
the  experimental  slack  (not  to  be  confused  with  bundle  slack)  was  very  touchy. 
When  toggling  the  grips  at  the  control  console  the  load  read  out  had  to  be 
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acutely  monitored  for  the  bundle  loading  must  be  kept  to  a  minimum  prior  to 
testing. 

During  testing  a  real  time  graphical  display  of  the  load  vs  displacement 
was  used  for  monitoring  the  results,  with  this  experimental  irregularities  in  the 
several  of  the  samples  was  detected.  The  bundles  in  the  very  short  gauge 
lengths  showed  irregular  load  curves.  Keying  into  this  after  several  samples 
helped  to  identify  a  problem  with  the  sample  preparation. 

Friction  was  expected  to  affect  the  long  gauge  length  (250  and  500  mm) 
and  not  affect  the  shorter  gauge  lengths  (25  and  50  mm).  To  test  the  affect 
each  gauge  length  tested  a  pair  was  treated  with  oil  and  a  pair  was  left 
untreated  (dry).  Tests  were  run  on  samples  of  500  mm,  250  mm,  50  mm, 
25  mm  and  5  mm. 

B.  DATA  REDUCTION 

After  tests  were  complete  the  data  files  were  decoded  into  ASCII  and 
stored  for  processing.  The  data  output  from  the  Instron  software  is  a  three 
column  file.  The  first  column  is  the  record  number,  the  second  is  the 
displacement  (inches)  and  the  third  is  the  load  (lbm).  This  is  then  processed 
through  a  program  to  convert  the  data  to  mm  and  kg  and  remove  compliance. 
There  are  two  data  files  output,  one  file  with  compliance  removed  and  the  other 
converted  raw  data  (compliance  not  removed)  [Ref.  5].  The  output  from  the 
conversion  program  is  then  processed  through  an  analysis  program  listed  in 
Appendix  C.  This  outputs  the  ECDF,  the  bundle  Modulus  E  and  the  Slack 
region  ECDF.  The  ECDF  is  then  processed  through  a  Weibull  Maximum 
Likelihood  Estimator  program  to  define  the  shape  parameter  alpha  and  the 
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location  parameter  beta.  ECDF  is  also  processed  through  a  linearizing  program 
which  outputs  the  Weibull  probability  plot  data  file. 


V.  RESULTS 


A-  TEST  OBSERVATIONS 

During  testing  several  of  the  specimens  were  invalid  due  to  experimental 
difficulties.  Some  of  the  samples  had  maximum  strains  noticeably  beyond 
expected  limits.  One  bundle  was  observed  to  have  a  much  lower  modulus  than 
what  was  expected  from  equation  C.5.  Dissection  of  the  samples  with  suspect 
loading  curves  revealed  a  lack  of  adhesive  wetting  resulting  in  a  large 
percentage  of  the  fibers  slipping  from  the  tab.  The  problem  appeared  most 
frequently  in  the  small  gauge  length  (GL)  Samples.  It  is  suspected  that  the 
adhesive  which  couples  the  bundle  and  tab  was  not  wetting  the  internal  fibers 
of  the  bundle,  with  a  resulting  loss  of  ability  to  carry  the  higher  loads  of  the 
small  GL.  This  phenomenon  manifested  itself  two  ways.  If  there  was  partial 
slipping  after  initiation  of  the  failure  region,  a  strain  higher  than  the  actual 
maximum  strain  was  observed.  If  total  slipping  occurred  in  a  large  group  of 
fibers  the  modulus  decreased.  After  closer  inspection  of  all  the  samples  it  was 
seen  that  many  of  the  specimens  had  a  small  number  fibers  slipping  in  the  tabs 
but  the  percentage  was  small  enough  not  to  effect  the  results. 

To  avoid  this  problem  the  remaining  tests  were  conducted  in  samples  at  or 
beyond  25  mm.  The  problem  was  less  acute  in  the  long  GL  because  of  the 
lower  loads  these  long  GL’s  experience  as  a  consequence  of  size  effect. 

Two  types  of  test  were  carried  out,  with  oil  and  dry.  This  data  was  then 
reduced  by  categories  of  compliance  removed  and  not  removed.  The  compliance 
was  estimated  by  testing  samples  of  zero  GL  [Ref.  5].  The  resulting  curves 
were  fitted  and  reduced  to  am  analytical  equation,  displacement  as  a  function  of 


load.  This  equation  can  then  be  used  to  subtract  the  compliance  from  the 
empirical  data.  The  zero  GL  samples  used  the  same  coupling  adhesive  and 
possibly  suffered  from  partial  slipping.  This  causes  the  compliance  curve 
coefficients  to  be  over  estimated  and  as  a  result  displacements  removed  for 
compliance  were  too  large.  When  the  compliance  was  removed  from  the 
experimental  data  displacements  in  the  lower  range  showed  negative  values 
indicating  over  estimations. 

Since  the  determination  of  compliance  is  slightly  in  error  it  is  not  used  in 
the  results  discussion,  though  data  reduction  calculations  were  completed  to 
indicate  possible  trends  and  affects.  Table  II  shows  the  data  reduction  with 
removal.  Table  III  is  without  compliance  removal  and  Table  IV,  for  comparison, 
is  the  result  from  single  fiber  tests  [Ref.  61. 


Table  II.  EXPERIMENTAL  DATA  RESULTS,  COMPLIANCE  REMOVED. 


Compliance  Removed 

! 

j 

Test  # 

Gage 

Modulus 

alpha 

Beta 

• 

dry/oil 

length 

(kg/mm) 

(mm) 

(gm/e) 

(mm/mm) 

100901 

25 

109.76 

915 

3.0 

0.0157 

dry 

090901 

50 

64.09 

1068 

3.214 

0.0135 

* 

dry 

100908 

250 

11.0 

917 

3.18 

0.00972 

dry 

100904 

500 

5.433 

905 

3.035 

0.0076 

dry 

9! 

090909 

50 

54.7 

912 

4.465 

0.0145 

oil 

100907 

250 

11.18 

932 

3.964 

0.0101 

oil 

100905 

500 

5.755 

959 

3.66 

0.0085 

oil 

w 
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B.  DATA  COMPARISON 


Generally  the  results  indicate  only  small  deviations  in  the  alpha  parameter 
(Table  III)  for  all  GL.  Between  a  oiled  sample  of  the  same  GL  only  the  500 
mm  GL  had  any  noticeable  change  in  the  shape  parameter,  which  had  been 
expected.  The  location  parameter  was  effected  most  by  the  oil  treatment,  beta’s 
increased  by  10%  in  both  the  50  and  500  mm  GL.  In  250  mm  GL  the  alpha 
and  the  beta  parameters  seem  to  be  relatively  stable,  meaning  the  oil  treatment 
had  little  effect. 

It  was  expected  that  friction  would  affect  the  longer  GL’s  and  have  little 
or  no  affect  on  the  shorter  samples.  The  500  mm  test  did  show  a  considerable 
difference  when  oil  was  applied  (Figure  5). 

The  friction  was  expected  to  become  influential  only  as  the  percentage  of 
failed  fibers  increased,  but  the  dry  bundle  shows  a  decreased  load  early,  which 
can  be  seen  in  the  ECDF  plot  (Figure  6).  In  fact  the  oil  and  dry  sample  data 
merge  as  the  percentage  of  failures  increases.  The  reduced  friction  made  a 
difference  in  the  resulting  Weibull  parameters,  alpha  increased  from  3.085  to 
3.8.  From  single  fiber  testing  alpha  (Table  IV)  for  AS-4  is  approximately  4.11. 
Even  with  oil  the  bundle  is  affected,  as  seen  by  the  sudden  drops  in  the  bundle 
load  (Figure  5)  or  the  sharp  increase  in  the  percentage  of  failures  (Figure  6)  at 
discrete  points. 

Without  friction  the  failure  distributions  are  expected  to  be  smooth. 
Decreasing  the  GL  it  is  anticipated  that  the  effects  of  friction  diminishes.  Tests 
at  250  mm  the  loading  curve  is  seen  to  smooth  out  (Figure  7),  providing  a 
slightly  more  continuous  plot.  Friction  had  less  effect  on  this  GL  though  a 
slight  decrease  is  noticed  when  the  bundle  is  dry. 
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Table  III.  DATA  REDUCTION,  COMPLIANCE  NOT  REMOVED. 


Compliance  Not  Removed 


Test  # 
dry/oil 

Gage 

length 

(mm) 

Modulus 

(kg/mm) 

(gm/e) 

Alpha 

Beta 

(mm/mm) 

Beta 
sized  to 
50  mm 

090901 

dry 

50 

45.6 

761 

4.75 

0.01615 

0.01615 

100908 

dry 

250 

9.85 

820 

4.07 

0.0104 

0.0154 

100904 

dry 

500 

5.2 

867 

3.085 

0.00781 

0.0137 

090909 

oil 

50 

43.5 

725 

4.89 

0.0179 

0.0179 

100907 

oil 

250 

10.6 

8  33 

3.95 

0.01041 

0.0154 

100905 

oil 

500 

5.5 

917 

3.8 

0.00868 

0.0152 

SUMMARY 

a 

Std 

13 

Std. 

E 

Std 

OIL 

4.21 

0.59 

0.0162 

0.0015 

825 

96 

DRY 

3.96 

0.84 

0.0151 

0.00125 

816 

53 

Merged 

4.09 

0.66 

0.0157 

0.00137 

821 

75 

Table  IV.  RESULTS  FROM  SINGLE  FIBER  TEST  AT  50  mm  GL. 

Fiber  Load  Weibull  Parameters 


Spool  # 

Alpha 

Beta 

008 

4.28 

0.0175 

019 

3.94 

0.0172 

average 

4.11 

0.0173 

Modulus  =  900  gmJe 
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Figure  5.  500  mm  gage  length  load  curve  comparison. 
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The  ECDF  plot  (Figure  8)  indicates  friction  had  more  influence  than  was 
revealed  in  the  load  curve.  The  upper  tails  again  merge,  indicating  either  the 
oil  no  longer  decreases  the  friction  or  friction  no  longer  affects  the  remaining 
fibers  in  this  region  of  the  ECDF.  A  slight  change  can  be  seen  between  the 
two  samples  in  the  lower  and  mid  region  but  this  could  be  variations  between 
bundles  and  not  whether  the  specimen  was  treated  with  oil,  but  it  should  not 
be  dismissed  because  it  is  consistent  with  the  other  samples. 

A  50  mm  GL  test  (Figure  9)  showed  results  which  had  been  anticipated  in 
the  longer  lengths.  Little  change  is  observed  as  few  fibers  fail  but  the  affects  of 
friction  become  prominent  as  the  percentage  increases. 

The  ECDF  (Figure  10)  indicates  the  effect  of  oil  more  clearly.  The  entire 
range  of  failure  is  shifted  not  only  at  the  upper  tail  but  the  lower  tail  as  well. 

Determining  how  well  the  underlying  Weibull  distribution  is  extracted  from 
a  bundle  test  and  what  influence  the  oil  treatment  has  can  be  observed  in 
Weibull  plots.  For  the  500  mm  sample  (Figure  11  &  12)  friction  has  definitely 
affected  the  distribution,  but  it  is  also  observed  that  oil  helped  to  reduce  some 
of  the  deviation  especially  in  the  lower  tail. 

As  the  gauge  length  is  reduced  it  becomes  apparent  from  the  Weibull  plots 
(Figures  13,14,15  &  16)  for  the  250  mm  GL  and  the  50  mm  GL  that  friction  is 
not  necessarily  negligible,  for  the  distribution  is  affected.  This  is  clearly 
demonstrated  in  the  50  mm  GL  treated  with  oil  (Figure  15)  which  shows  the 
characteristic  Weibull  distribution  distinctly  as  apposed  to  the  dry  sample. 
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Figure  7.  250  mm  gage  length  load  curve  comparison. 
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ECDF  for  50  mm  comparison. 
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Figure  11.  500  mm  GL  Weibull  plot,  oil  treated. 


2.50 


29 


30 


2.50  -2.00 


31 


2.50  -2.00  -1.50  -1.00  -0.50  0.00  0.50 

L0G(X/BETA) 


2.00 


33 


2.50  -2.00 


VI.  SUMMARY  AND  RECOMMENDATIONS 


A.  SUMMARY 

The  affect  of  friction  is  evidenced  in  all  gauge  lengths  and  is  not  always 
limited  to  the  upper  tail.  Lubrication  of  the  bundle  reduces  the  effect  of  friction 
and  should  be  used  in  all  bundle  failure  testing. 

The  compliance  removal  was  not  used  in  the  data  reduction  due  to  the 
error  in  the  compliance  curve  caused  by  adhesive  coupling  failure  in  the  tabs. 
The  compliance  ended  up  being  over  estimated.  This  a  side  the  trend  of 
compliance  removal  on  the  data  reduction  (Table  II)  indicated  a  general 
reduction  in  the  parameters  (a,G).  Estimating  the  compliance  to  be  much  less 
and  extrapolating  this  thought,  it  can  be  assumed  the  value  for  the  parameters 
presently  determined  would  be  slightly  less.  In  a  comparison  of  the  single 
fibers  parameters  with  the  bundle  parameters  (Table  III),  the  bundle 
parameters  were  barely  higher.  Given  the  affect  of  compliance  it  can  be 
concluded  that  once  removed  the  bundle  parameters  will  be  indistinguishable 
from  the  single  fiber  testing  and  in  fact  given  the  quantity  of  the  statistical 
base  much  more  reliable. 

It  was  shown  that  friction  effected  all  the  gauge  lengths  but  the  longer 
gauge  length  (500  mm)  experienced  the  most  disturbance.  The  50  mm  gauge 
lengths  when  treated  with  oil  produced  a  very  smooth  distribution  in  the 
Weibull  plot.  The  dry  50  mm  sample  produced  a  poor  Weibull  plot.  The  50 
mm  gauge  length  produce  good  results  through  out  the  entire  range.  The  500 
mm  gauge  length  produced  good  results  in  the  Weibull  plot  at  the  lower  tail 
(Figure  7).  Since  the  lower  tail  of  the  large  gauge  length  are  statistics  of  much 
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higher  reliability  than  the  50  mm  the  linear  statistics  of  the  500  mm  can  be 
pieced  into  the  50  mm  statistical  plot  there  by  extending  the  range  of  the 
variable  (strain).  The  upper  tail  of  the  50  mm  gauge  length  is  not  completely 
linear  in  the  Weibull  but  this  can  be  removed  and  statistics  from  a  5  mm  test 
can  replace  it.  This  piece  meal  work  on  the  Weibull  plot  can  be  done  because 
of  size  effect.  The  slope  alpha  does  not  change  but  beta  decreases  in  the 
Weibull  as  the  gauge  length  increases.  With  this  procedure  friction  and  other 
contaminates  can  be  removed  and  the  underlying  parameters  accurately 
predicted. 

Through  simulation  a  normal  or  uniform  distribution  of  slack  was 
determined  to  have  very  little  effect  on  the  distribution  if  the  dO  displacement 
is  subtracted. 

B.  RECOMMENDATIONS 

To  ensure  better  results  several  areas  of  experimentation  can  be  improved 
on. 

-  The  handling  of  the  specimens  should  be  reduced,  a  rigid  carrying 
mechanism  which  can  transport  the  bundle  from  assembly  to  storage 
to  testing  without  any  human  handling  or  disturbing  the  fibers  is 
necessary. 

-  The  manual  mechanical  grips  should  be  replaced  with  a  pneumatic  device. 

-  Oil  treatment  should  be  used  on  all  samples. 

-  Improve  the  tab  coupling. 

The  slack  distribution  requires  more  investigation  than  was  addressed  here. 
Simulation  was  a  key  element  which  helped  to  determine  it’s  effect  on  the 
distribution.  An  analysis  algorithm  is  provided  in  this  thesis  which  will  help  to 
determine  the  slack  ECDF.  To  further  study  the  question  of  slack  several 
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bundle  samples  (large  gauge  length)  dedicated  solely  for  the  determination  of 
slack  should  be  tested.  A  much  larger  data  base  is  required  to  use  the  methods 
discussed  in  Appendix  C.  By  decreasing  the  cross-head  rate  and  increasing  the 
acquisition  rate  a  large  data  file  can  be  attained.  Once  the  distribution  is  fitted 
to  a  model,  simulation  can  determine  what  affect  it  has  on  the  bundle  strength 
distribution. 


vn.  CONCLUSIONS 


The  results  achieved  in  this  thesis  conclude  convincingly  that  the 
distribution  of  the  fiber  can  be  retrieved  through  bundle  failure  tests.  The 
parameters  extracted  from  the  testing  compared  favorably  with  the  single  fiber 
testing.  As  a  result,  enormous  characterization  efficiency  is  realized;  in  fact 
with  this  technique,  statistics  to  10*  are  realized. 
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APPENDIX  A 


SIZE  EFFECT 


The  determination  of  the  scale  parameter  beta  of  a  gage  length,  which  is  too 
short  to  preclude  practical  experimental  implementation,  can  be  accomplished  by 
size  effect  of  the  weibull  distribution  function. 

Starting  with  the  Weibull  CDF  for  a  fiber  of  defined  gage  length  1,: 


F,(x)  =  1  -  expKx/JJ.r') 


R,(x)  =  1  -  F,(x) 


For  a  second  fiber  of  gage  length  12  with, 


1,  <  lj  and  L/l,  =  m  ;  m  =  integer  >  0 


and  assuming  the  reliability  R,  is  constant. 


R,  =  (R,)(R,XR3)....(RJ  =  R, 


substituting  A.2  for  R, 


R,  =  exp{-(x/B,nr 
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and  the  total  reliability 


Ri  =  expMx/flJ-} 

equating  equation  A.4  and  A.  5 

expWx/B.r)  =  expKx/B,nm 
taking  the  natural  log  A.6  reduces  to 
(x/Q„r  =  m  (x/fl.r* 
again  taking  the  natural  log  of  A.7 
aJnCx/flJ  -  a,ln{m(x/fl,)l  =  0 
since  a,,  =  a,  =  a  is  not  =  0 

BA  =  m*'a 
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APPENDIX  B 


EXPERIMENTAL  PROCEDURES 

A.  EXPERIMENTAL  SET  UP 

The  Test  apparatus  consisted  of  an  INSTRON  Model  4206  material  tester, 
with  a  50.0  kg  load  cell,  connected  to  a  4200  Series  Expanded  Control  Console. 
The  Console  was  connected  to  a  IBM  PC/AT  through  com  port  1  via  a  IEEE 
connection  and  converter.  The  data  acquisition  and  material  tester  control  was 
commanded  by  Instron  Series  IX  Automated  Material  Testing  System  series 
4.01C  software. 

The  controlling  software  has  a  5  page  menu  which  requires  setting  prior  to 
testing.  Many  of  these  settings  are  for  Instron  data  reduction  processing  which 
are  not  used  in  this  thesis  but  required  values  to  operate.  The  settings  crucial 
to  testing  and  proper  acquisition  are  cross-head  speed,  data  acquisition  rate  and 
gage  length.  Safety  features  include  a  maximum  load  and  displacement  setting 
to  insure  the  load  cell  and  testing  model  are  not  damaged. 

B.  INSTRON  CALIBRATION  AND  TESTING 

The  Instron  4206  and  control  unit  required  a  warm  up  time  of  one  hour 
prior  to  operation  to  allow  for  stabilization.  The  load  cell  has  two  methods  for 
calibration,  an  electronic  calibration  and  a  mechanical  calibration.  A  mechanical 
calibration  was  performed  for  all  tests. 

1.  Mechanical  Calibration 

Any  inputs  refer  to  the  4200  Series  Expanded  Control  console. 

-  Main  power  on/off/on  -  wait  for  diagnostic 
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-  Press  Load  Balance  /  Enter 

-  Hand  prescribe  a  5.0  kg  weight 

-  Press  Load  Cal 

-  Enter  the  weight  /  Enter 

-  Remove  the  weight 

-  Enter  Load  Balance  /  Enter 

-  Re-hang  the  weight  to  verify  correct  calibration 

-  Repeat  if  required 

2.  Loading  the  Sample 

The  procedures  for  sample  placement  into  the  test  apparatus  are  as 
follows. 

-  The  Instron  specimen  grips  are  separated  enough  to 
allow  the  sample  to  hang  from  the  top  grip  without 
touching  bottom  grip.  The  tab  of  the  specimen  is 
firmly  clasped  with  tweezers  and  gently  lifted  off 
resting  bench  allowing  to  hang  free. 

-  The  samples  top  tab  is  placed  in  the  top  grip, 
careful  to  align  the  bundle  with  the  center,  which 
has  been  measured  and  marked.  After  alignment  the 
grip  is  tightened  making  sure  not  to  twist  or  damage 
the  specimen. 

-  A  load  balance  is  performed  at  the  control  panel  to 
remove  the  weight  of  the  specimen. 

-  The  cross-head  is  then  toggle  down  so  as  to  place 
the  bottom  tab  in  the  bottom  grip  then  tighten. 

Note:  after  tightening  the  grips  the  specimen  develops 
noticeable  slack  and  the  load  cell  measures  a  negative 
load.  The  specimen  fibers  are  slightly  compressed  and 
need  to  be  straightened. 

-  The  cross-head  is  toggle  up  until  a  load  of  0.1  kg 
is  indicated  in  the  load  readout  window  on  the 
control  panel.  The  load  is  then  brought  down  to 
zero  and  the  displacement  reference  zeroized.  The 
procedure  is  repeated  once  to  ensure  the 
displacement  returns  to  at  or  near  zero.  The 
specimen  is  ready  for  testing. 

3.  Testing 

The  test  is  initiated  at  the  IBM  PC/AT  through  the  Instron  software. 
After  the  specimen  test  parameters  are  set  and  loaded,  prior  to  test  initiation, 


the  IEEE  port  is  enabled.  The  software  then  signals  the  system  is  ready  to 
begin  testing. 

As  the  cross-head  displaces,  load  and  displacement  data  are  sent  to  the 
PC/AT  at  regular  intervals  of  displacement,  as  prescribed  by  the  initial  software 
settings. 

A  real  time  graphical  output  is  available  to  monitor  the  testing.  This 
allowed  for  a  quick  identification  of  bad  tests.  After  the  test  was  complete,  the 
data  file  is  converted  from  Instron  System  code  to  ASCII  for  latter  data 
analysis. 

C.  SAMPLE  PREPARATION 

The  samples  were  made  from  a  Hercules  Magnamite  high  strength  graphite, 
type  AS-4,  spool  145.  The  bundle  has  3000  fibers  with  a  denier  of  .0057446 
gm/in.. 

1.  Procedures 

A  length  of  bundle  is  clamped  at  one  end  on  a  3  meter  aluminum  track 
and  weighted  on  the  other.  Half  Copper  tabs  approximately  2.54  cm  (width)  x 
2.54  cm  (length)  are  pre-positioned  on  the  track  at  a  specified  gage  length. 
Slots  built  into  the  track,  center  the  bundle.  An  adhesive  is  applied  at  the  tab, 
the  bundle  and  the  remaining  half  of  the  copper  tab  are  securely  clamped  on  to 
complete  the  tab.  The  bundle  is  severed  between  samples  on  the  track  and  the 
adhesive  is  allowed  to  cure. 

Once  dry  the  samples  are  manually  removed  from  the  track  and  placed 
on  the  sample  bench. 


■ 


i 


i 
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D.  CALCULATION  OF  CROSS-HEAD  SPEED 

The  cross-head  calculation  is  a  function  of  the  number  of  data  points 
required,  data  acquisition  rate,  load  cell  and  displacement  gage  tolerance. 

The  time  to  bundle  failure  was  selected  to  remain  constant.  This  will 
provide  a  constant  strain  rate  on  the  different  gage  lengths  that  were  tested. 

The  relationship  between  the  cross-head,  acquisition  rate  and  recorded  range 
is  given  in  equation  B.l. 

XH  =  (RRXARVDPTS  (B.l) 

where:  XH  =  cross-head  speed 

RR  =  recorded  range 
AR  =  acquisition  rate 
DPTS  =  number  of  data  points  recorded 

The  sampling  interval  is 

deltaD  =  RR  /  DPTS  =  XH  /  AR  (B.2) 

In  order  not  to  sample  beyond  the  abilities  of  the  acquisition  system  the 
following  rule  must  be  followed 

deltaD  >  Recording  Tolerance  (B.3) 

The  Instron  output  tolerance  for  displacement  was  0.00254  mm,  therefore 
deltaD  >  0.00254  mm.  Using  this  as  a  guideline  and  the  constant  testing  time 
chosen  of  5  min..  Table  IB  lists  the  Instron  test  setting  calculations. 
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Table  IB.  TEST  SETTINGS. 


Gage 

length 

(mm) 

Max 

displ. 

(mm) 

Time 

(min) 

Cross- 

head 

(mm/s) 

Acq. 

rate 

(pts/s) 

500 

25 

5 

8.33xl02 

6.6 

250 

12.5 

5 

4.17x10  s 

3.3 

50 

2.5 

5 

8.33x10  s 

3.3 

25 

1.25 

5 

4.17x10  s 

1.67 
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APPENDIX  C 
DATA  ANALYSIS 

The  analysis  of  the  bundle  failure  load  vs  displacement  curve  can  be  divided 
into  three  distinct  regions  (see  Figure  1C): 

-  Slack  region. 

-  Linear  region. 

-  Non-linear  region. 

Slack  is  created  during  specimen  production  and  placement  in  the  Material 
testing  machine.  As  the  bundle  is  measured,  cut  and  tabs  are  placed  on  each 


I 

I 

) 

end  of  the  sample  some  of  the  individual  fibers  become  loose  and  effectively 
i  longer  than  other  fibers.  This  produces  an  uneven  loading  across  the  fibers  as 

the  cross-head  of  the  tensile  testing  machine  strains  the  bundle.  The  ilh  fiber  is 


LOAD  CURVE 

REGIONS 


'igure  1C.  Regions  of  the  load  curve. 


t 
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not  loaded  until  the  displacement  of  the  cross-head  is  equal  to  the  respective 
slack  of  the  i*  fiber. 

When  the  cross-head  displaces  beyond  the  maximum  slack,  the  loading  curve 
becomes  linear  until  the  first  failure. 

The  amount  of  slack  varies  according  to  an  undefined  distribution.  The 
slack  distorts  the  failure  region  thereby  disturbing  the  underlying  sequence  of 
the  strength  distribution.  The  effect  slack  has  on  the  underlying  shape 
parameter  depends  on  the  slack  distribution  and  amount  of  slack  relative  to  the 
total  strain,  this  is  demonstrated  in  Appendix  D. 


A.  SLACK  DISTRIBUTION 

The  slack  in  a  bundle  is  dependent  on  the  variation  in  the  fiber  lengths. 


L*  =  1  +  delL,  (C.l) 

where:  Lj  =  length  of  fiber  j 

1  =  mean  gage  length 
delLj  =  slack  of  fiber  j 

Since  the  delL,  «  1  the  length  is  approximated  by 

1^  =  1  (C.2) 

The  load  in  the  slack  region  can  be  defined  by  Hooke’s  law. 

k-l 

Ps(dA)k  =  I  E  t(dA  -  delLj)/(l  +  dell^)]  (C.3) 

j=i 


k-l 

«  IE  d/1  ;  k=l,2,....,n 

j=i 

where:  Ps(  \  =  bundle  load  as  each  fiber  k  is  loaded 

E  =  fiber  modulus  (gm/strain) 
dA  =  apparent  displacement  =  d>  +  delL, 
d,  =  fiber  displacement 
n  =  number  of  fibers 


Figure  2C.  Histogram  representation  of  a  PDF  plot. 

Defining  the  loading  slack  in  terms  of  the  distribution  function  (fs)  of  Ps  at 
a  displacement  dA  as  shown  in  Figure  2C.  Ps  can  be  write  as 


Ps(d),  =  {fs(dA1XdA2-dA1)n)  E  U(d-dAI ){(d-dAI )/!} 


where:  d  =  any  displacement  such  that  d  >  dA) 

U(d)  =  unit  step  function  {=1  d  >  dA); 

{=0  d  <  dA1 

Summing  the  load  (eqn.  C.4)  over  the  range  of  dA  and  defining  the  bundle 
modulus  as  a  function  of  the  fiber  modulus 


E  =  nE/1  (C.5) 

Ps(d)  =  E  I  (fs(dA  jXdA  M  -  dA ,)]  U(d-dA ,)  (C.6) 
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differentiating  with  respect  to  d  for  a  continuous  function 

fd, 

8P/8(d)  =  E  fs(d)  d(d)  =  E  F(d.)  (C.7, 

Jd, 

F(d,)  is  the  CDF  of  the  slack.  Assuming  the  empirical  data  approximates 
the  continuous  data. 

P*  =  F  -  E  =  F  (C.8) 

The  derivative  of  the  empirical  data  can  be  numerically  solved  with  two 
methods,  a  numerical  differentiation  or  a  discretization. 

1.  Numerical  Differentiation 

Even  interval  Forward  and  Central  difference  methods  are  used  due  to 
the  nature  of  the  slack  loading  and  the  recording  system.  As  the  fibers  load 
the  slope  is  discontinuous  as  seen  in  Figure  3C.  Several  numerical  methods 
exist  which  will  produce  varying 
accuracy.  In  trying  to  recover  the 
underlying  distribution  of  the  slack 
the  sampling  rate  and  number  of 
fibers  must  be  considered  (i.e. ,  how 
many  data  points  can  be  recorded 
relative  to  the  number  of  fibers). 

If  the  sampling  rate  is  fast 

and  the  number  of  total  data  points  i 

recorded  in  the  slack  region  is  2-3 

times  that  of  the  fiber,  a  3  point  Figure  3C.  Slack  region  showing 

discrete  data  on  the  underlying  loading 
slope. 


forward  difference  method  may  yield  accurate  results  due  to  the  linearity 
between  discrete  data.  If  the  relative  sampling  rate  is  slow  the  3  Pt.  forward 
difference  may  skew  the  distribution.  To  prevent  this  a  central  difference 
method  is  safe. 

3  pt.  Central  Difference  method. 

E’  =  (-Em  +  Ei.ll/2h  (C.9.1) 

5  pt.  Central  Difference  method. 

EXE  ,2-8E  „-8E  j.rE  uaV2h  (C.9.2) 

where:  h  =  d  Al.,  -  dAl 

2.  Discretization 

The  data  in  the  slack  region  can  be  fit  to  a  curve  and  a  equation 
estimated.  This  equation  can  be  differentiated  analytically.  The  underlying 
data  itself  is  not  smooth  so  such  a  method  is  not  recommended. 

B.  LINEAR  REGION 

The  slope  of  the  linear  region  is  the  bundle  modulus.  To  calculate  this 
slope  a  least  squares  method  is  very  accurate,  due  in  part  to  the  natural 
linearity  of  the  data.  Before  this  method  can  be  employed  the  region  must  be 
specifically  defined.  The  point  where  the  non-linear  slack  region  ends  is 
identified  as  dl  and  the  end  of  the  linear  region  is  identified  as  d2  (the  point  of 
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the  first  failure).  This  is  graphically  displayed  in  Figure  AC.  Between  dl  and 
d2  the  least  squares  method  calculates  modulus  £. 

The  points  dl  and  d2  must  be  graphically  identified  and  manually  entered 
into  the  analysis  program.  To  do  this  a  plot  of  the  load  curve  must  be 
generated  before  the  analysis  can  be  conducted. 


LOAD  CURVE 


Figure  4C.  The  linear  region. 
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1.  Least  Squares  Method 

Due  to  the  noise,  the  discrete  data  i  will  deviate  slightly  from  the 
underlying  slope  £.  From  Hooke’s  law  the  expected  relationship  is 


Ei  =  -ED, 


(C.10) 


where:  D,  =  dh-dl 


(C.10.1) 


Figure  5C.  Least  Squares  Method. 


But  due  to  noise  the  load  £  deviates  by  an  amount  Pe,  this  is 
graphically  represented  in  Figure  5. 

Pej  =  ED,  -  P,  (C.ll) 

Defining  the  least  squares  sum  M 

M(E)  =  I  ((ED,  -  El/202)2  (C.12) 

where:  o  =  standard  deviation 

Differentiating  M  with  respect  to  E  and  setting  it  equal  zero  defines  the 
minimum  least  squares. 

dM/dE  =  0 
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substituting  in  C.10.1  this  reduces  to 

E=  lE.tlLrdD/KiL-dl)1  (C.13) 

2.  Intercept  of  the  Abscissas 

With  the  Modulus  E  the  intercept  of  the  abscissas  can  be  calculated. 

dO  =  dl  -  E(diy£  (C.14) 

C.  FAILURE  REGION 

This  region  is  the  key  to  retrieving  the  strength  distribution.  Using  the 
slope  of  the  bundle  the  failure  statistics  can  be  transferred  to  a  ECDF  space. 
From  the  ECDF  the  parameters  of  the  distribution  can  be  estimated. 

The  ECDF  is  a  plot  of  the  percentage  of  failures  versus  a  variable. 
Typically  this  variable  is  the  strength,  this  would  be  difficult  to  translate  due  to 
the  nature  of  the  bundle  test.  The  displacement  is  constantly  incremented, 
where  as  the  load  is  a  bundle  load  and  can  vary.  Because  of  this  the  ECDF  is 
plotted  versus  the  displacement.  The  displacement  is  proportional  to  the 
strength  as  defined  by  Hooke’s  law  and  once  the  fiber  modulus  is  identified  the 
displacement  can  be  translated  to  load. 

When  a  fiber  fails,  the  bundle  modulus  E  changes  by  one  fiber  strength. 
From  Hooke’s  law  it  is  expected  that  the  failures  would  occur  along  the 
modulus  Ej.  Ej  being  the  modulus  of  the  unfailed  bundle  which  ideally  is  a 
multiple  of  the  single  fiber  modulus  proportional  to  the  number  of  remaining 
fibers.  Unfortunately  this  is  not  the  case,  as  the  bundle  loads,  friction  and 
slack  change  the  distribution  causing  the  actual  loading  path  to  deviate  from  E; 
and  follow  along  £;• 
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To  show  this  in  more  detail  consider  first  the  ideal  underlying  loading  (no 
slack,  friction  or  noise)  of  a  bundle  under  constant  strain.  At  a  break  in  one 
fiber  j,  the  load  will  drop  until  intercepting  the  slope  of  shown  in  Figure  6C. 
This  new  slope  emanates  radially  out  from  the  intercept  of  the  abscissas,  (note: 
without  slack  the  intercept  should  be  the  origin).  If  the  data  acquisition  rate 
where  high  enough  this  failure  could  be  clearly  identified.  Now  add  slack,  this 
extends  the  failure 
displacement  of  some 
fibers  and  the  intercept 
of  the  abscissas  by  E, 
after  the  initial  failure 
to  some  point  between 
the  origin  and  the  dO 
point.  The  intercept 
will  shift  a  relatively 
small  amount  but  will 
digress  toward  the 
origin.  If  friction  is 
present  the  broken  fibers  will  cling  to  the  unbroken  fibers  and  cause  additional 
strain.  Add  noise  to  the  acquisition  system  and  the  discrete  data  fluctuates. 
These  factors  combine  to  force  the  actual  recorded  data  to  deviate  from  the 
underlying  bundle  loading  path. 

A  bundle  with  n  number  of  fibers  will  have  n  failures.  If  the  data 
acquisition  system  can  not  record  this  many  points  in  the  failure  displacement 
range  than  each  statistic  can  not  be  recorded.  With  slack,  friction  and  noise 


Figure  6C.  Failure  loading  path. 
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even  with  a  fast  enough  data  acquisition  system  the  failures  might  not  be 
distinguishable. 

The  percentage  of  failures  in  the  bundle  test  can  be  defined  by  the  ratio  of 
the  modulus  of  no  failure  to  the  modulus  of  failures.  The  reliability  is  defined 
by 

&  =  E/E  =  E(dAlVlE(dAi-dO)l  (C.15) 

The  percent  failure  is  1  minus  reliability,  the  ECDF  is  defined  by 

EOU  =  1  -  fi,  (C.16) 

This  relationship  is  graphically  depicted  in  Figure  6C.  The  ECDF  begins  at 


Load 

P 


ECDF 
F  (%) 


Figure  7C.  Transfer  from  loading  curve  to  ECDF  space. 
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62,  so  the  manual  selection  of  this  point  must  be  done  with  care. 

Assuming  all  conditions  idea!  the  ECDF  would  fall  out  as 

F  =  j/n  (C.17) 

But  because  conditions  are  not  ideal  this  does  not  occur.  What  is  important  is 
the  shape  of  the  ECDF  and  whether  it  is  close  to  the  underlying  CDF.  The 
ability  to  distinguish  each  failure  is  unimportant.  Another  point  to  note  is  the 
lower  tail  needs  to  be  clearly  defined.  Due  to  the  relatively  low  failure  rate  in 
the  lower  tail  and  a  constant  recording  rate  this  area  of  the  curve  has  a  high 
number  of  discrete  data  points  per  failure,  so  in  essence  the  testing  method 
takes  care  of  itself. 

1.  Upper  Bound 

The  data  is  recorded  on  even  intervals  of  displacement  not  just  at  the 
failure  points,  this  is  graphically  described  in  Figure  6.  The  discrete  data  most 
closely  resembling  the  failure  point  is  on  the  top  edge  of  the  loading  curve.  This 
is  referred  to  as  the  Upper  bound.  When  this  is  translated  into  the  ECDF 
space  the  Upper  bound  is  on  the  far 
edge  of  the  plot  as  seen  in  Figure  8C. 

To  retrieve  the  underlying  CDF  from 
the  ECDF  the  Upper  bound  (Fu) 
must  be  extracted. 

Attempting  to  extract  the 
Upper  bound  from  a  bundle  of  n 
fibers  the  most  failure  points  that 
could  be  optimistically  retrieved  is  n. 

With  a  large  n  this  is  not  possible 
due  to  the  acquisition  rate,  no 


Figure  8C.  Close  view  of  the  ECDF. 
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matter,  the  shape  and  location  is  important  and  does  not  require  a  large 
number  of  points. 

Using  the  empirical  ranking  as  a  reference  point 

F  =  k^m  ;  k  =  1,2, . m  (C.18) 

where  m  is  the  number  of  points  desired  in  the  ECDF.  F  is  used  to  define 
search  bands  of  data  where  the  upper  and  lower  band  is  defined  by  1/n.  The 
value  of  E,  which  has  the  largest  displacement  (d^)  is  assigned  to  Ell*  When 
no  discrete  data  is  in  the  search  band  F,  Euk  is  assigned  the  value  of  F  at  a 
displacement  dA  k 

2.  Maximum  Likelihood  Estimator. 

With  the  Upper  bound  extracted  this  data  can  be  used  to  estimate  the 
parameters  of  the  distribution.  The  model  selected  is  a  Two  Parameter  Weibull 
Distribution.  The  estimator  chosen  is  the  Maximum  Likelihood  Estimator 
(MLE)  [Ref.  #:p.  103]  developed  by  R.  A.  Fisher.  T'  ,  method  finds  the 
maximum  of  the  likelihood  of  the  sample  in  the  range  of  the  variable  dll  and  is 
a  function  of  distribution  parameters.  For  the  Weibull  model  the  method  solves 
for  the  two  parameters  alpha  and  beta.  The  parameter  beta  is  solved  explicitly 
[Ref.  4]  in  terms  of  alpha.  Alpha  can  not  be  solved  explicitly  and  an  iteration 
must  be  used. 

13  =  [  17m  I  du,2  ],<2  sum  i=l  to  m  (C.19) 


a  =  {  I  dll,2  In(du,)AT  dll2)  *  £  IrsCdii.)) 1 


(C.20) 


An  Atkinson  method  of  iteration  was  chosen,  initially  alpha  is  assumed.  From 
this  successive  alpha’s  are  calculated  using  equation  C.20. 

a’.  =  a.  -  (a„  -  a,)2/(a„  +  2a,  +  ctj  (C.21) 

where:  a,  =  initial  start  value 

a,  =  f„(a*) 
a?  =  fa(a,) 

This  equation  is  iterated  for  a’0  until  a’0  =  a„. 

1.  Linearized  ECDF 

The  data  is  linearized  for  a  Weibull  probability  plot  to  visually  check 
the  conformity  to  this  model.  To  linearize  the  ECDF  start  with  the  Weibull 
equation 

F  =  1  -  exp(-P/G)a 

subtracting  1,  multiplying  by  -1  and  taking  the  natural  log. 

ln(l-F)  =  (-P/J3)0 
again  taking  the  natural  log 

ln[-ln(l-F)]  =  aln(P/J3)  (C.22) 

C. 22  is  a  linear  line  with  a  slope  of  alpha. 

D.  DATA  OUTPUT 

As  the  analysis  program  calculates  certain  data  file  are  output,  these  files 

are: 

-  EXPER.OUT,  the  converted  experimental  output  displacement,  strain 
and  load. 

-  SLACK.OUT,  the  slack  region. 

-  EXMS.OUT,  the  load  curve  minus  slack  region. 
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-  UBFECDF.OUT,  the  Upper  bound  ECDF. 

-  FECDF.OUT,  the  ECDF. 

-  SECDF.OUT,  slack  ECDF. 
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APPENDIX  D 


A.  SIMULATION 

The  simulation  of  a  bundle  test  has  two  purposes.  First  by  analyzing  in 
depth  the  mechanisms  which  govern  the  failure  of  the  fibers  in  a  bundle,  a 
better  understanding  of  how  to  extract  the  characteristic  properties  is  attained. 
The  properties  desired  are  the  distribution  parameters. 

Secondly  a  simulation  gives  us  the  ability  to  know  what  the  underlying 
parameters  are  and  thereby  proof  the  analysis  data  reduction  algorithm.  Also 
the  variables  which  influence  the  shape  of  the  failure  distribution  can  be  varied 
in  simulation  and  the  degree  of  induced  affect  observed. 

The  bundle  loading  curve  can  be  divided  into  three  regions; 

-  Slack 

-  Linear 

*  Failure 

Initially  the  simulation  requires  input  of  the  variable  estimates,  this  is 
done  through  an  input  file. 

1.  Initial  Inputs 

The  underlying  parameters  define  the  strength  distribution  for  a  gage 
length.  The  distribution  chosen  is  the  Weibull  model.  The  parameters  are 
alpha  and  beta,  the  number  of  fibers  per  bundle  is  n.  The  fiber  modulus  (E)  is 
in  gm/mms,  from  Hooke’s  law 

o  =  E  e  (D.l) 

where:  o  =  stress 

e  =  strain 
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Substituting  the  definition  of  stress  and  strain  for  fiber  i 


P/A  =  E  8/Li  (D.2) 

where:  P  =  load 

A  =  area  =  constant 
L  =  length  of  fiber 
6  =  change  in  length  L 

rearranging  D.2  and  letting  Ef  =  EA  (gm/e) 

P,  =  Ef  5/L,  (D.3) 

2.  Slack  and  Strength  Distributions 

The  slack  is  caused  by  the  variation  in  fiber  length  L.  If  L  were 
constant  no  slack  would  be  evident.  The  slack  is  expected  to  have  a  definable 
distribution,  but  this  is  unknown  and  could  vary  from  bundle  to  bundle 
therefore  several  models  are  simulated.  Two  models  selected  where 

-  Two  parameter  Weibull 

-  Uniform. 

Although  physical  consideration  favors  the  slack  model  to  be  normal, 
the  Weibull  model  was  chosen  for  it’s  analytical  convenience.  The  Weibull  is 
very  flexible  and  can  characterize  a  wide  range  of  shapes.  For  an  alpha  of  3.5 
the  distribution  is  approximately  normal,  other  than  3.5  the  shape  of 
distribution  is  skewed  negative  or  positive.  Beta  identifies  the  central  tendency 
of  the  distribution,  its  relatable  to  the  mean.  From  equation  2.1  the  Weibull 
CDF  is  defined  as 

F„(x;a,(3)  =  1  -  exp(-(x/13)“ 
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rearranging  in  terms  of  the  variable  x,,  which  can  represent  slack  (delL,)  or 
strength  P(. 

x,  =  exp(  [ln(-ln(l-Fw)  +  cdn(B)]  /  a)  (D.4) 

0  <  F.  <  1 

The  uniform  model  is  a  simple  two  variable  distribution  with  a  lower 
(u,)  and  upper  (u2)  range.  Since  the  slack  is  initiated  at  zero  the  lower  range  is 
zero. 

Fu(delL;u,,u2)  =  (delL  -  u,)/(u2  -  u,)  (D.5) 

0  <  Fu  <  1 

solving  for  delL,  from  D.5  and  setting  u,  =  0 

delL,  =  u2  Fu  (D.6) 

3.  Data  Generation  and  Continuity 

A  uniform  random  number  generator  creates  an  array  of  numbers 
between  0  and  1  and  stores  it  until  recalled. 

Recalling  equation  D.3  and  substituting  in  C.l  and  solving  for  the 
displacement  (8,). 

8,  =  P,  (1  +  delL,)  /  Ef  (D.7) 

The  bundle  is  put  under  a  constant  strain  rate,  so  the  discrete 
variable  will  be  displacement.  As  the  bundle  is  strained  the  filaments  will  fail 
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in  order  of  the  underlying  strength  or  strain  (neglecting  slack  and  friction).  For 
a  constant  filament  length  the  strain  is  proportional  to  the  displacement  and 
the  fiber  will  fail  at  some  underlying  displacement.  If  slack  is  present  in  the 
bundle  each  fiber  will  not  fail  at  the  underlying  displacement  because  the 
displacement  is  no  longer  proportional  to  strain  but  a  function  of  both  slack  and 
stain. 

The  displacement  observed  by  the  Instron  material  testing  machine  is 
not  5.  The  machine  records  an  apparent  displacement  (dA) 

dA  j  =  5,  +  delL,  (D.8) 

This  will  affect  the  fiber  failure  order,  every  other  variable  is  then 
keyed  to  the  order  (r^)  of  dAj. 

4.  Slack  Loading 

As  the  bundle  is  strained,  the  fibers  are  loaded  by  a  progressive 
summation  of  each  fiber  in  order  of  their  slack  distribution.  The  loading 
function  in  the  slack  region  can  be  defined  as  a  function  of  the  displacement  (8.) 
each  fiber  experiences  after  slack  is  released. 

j  i 

P.j=  £Ef5./U  ;  j-  1,2 . ,n  (D.9) 

i=l 

j-1 

P.j  =  E  Ef  [(delLy,  •  delLo  ,)/(l  +  delLo ,)]  (D.  10) 

i  =  l 

where:  delL*,,  =  order  slack 
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This  relationship  can  be  represented  graphically  in  Figure  ID.  The  last  load 
calculated  from  equation  D.10  will  occur  at  delL> 


Ps3 

Pa 

P61 

delL  delL  delL0 

1  2  3 

Figure  ID.  Slack  region  load  simulation. 

4.  Failure  loading 

The  order  of  each  fiber  failure  is  determined  by  the  ordered  strength 
distribution  (assuming  no  slack).  With  slack  the  order  of  failure  will  deviate 
slightly.  The  failure  load  is  a  function  of  the  modulus  of  the  bundle  and  the 
displacement  of  each  fiber.  The  bundle  failure  load  (Pu)  can  be  defined  from 
equation  D.2  as 


Puj  =  X  Ef  6/L 

substituting  in  C.l  and  D.8 


PU)  =  X  E  [(dAj  -  delL, V(1  +  delL] 
H 


(D.ll) 


(D.12) 
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5.  Discrete  Data 

The  material  testing  machine  measures  and  transmits  data  on  even 
intervals  of  displacement.  To  simulate  this  the  defined  loading  points 
(P„dell<j)  and  (Pu,dA)  will  be  set  limits  for  a  linear  interpolation.  The  data  is 
interpolated  on  intervals  of  deldA  which  is  set  by  the  cross  head  speed. 
Defining  the  discrete  apparent  displacement 

£Lu!=  & ,  +  deldA  (D.13) 

The  interpolation  algorithm  to  define  the  discrete  data  from 
continuous  is  the  same  for  both  slack  and  linear  region. 

The  discrete  load  P  in  the  slack  region  is  a  function  of  P.,  the  ordered 
slack  (delL0)  and  the  apparent  discrete  displacement 

E  =  KP.,.,  -  P,  jlAdelLo,.,  -  delL0j)]dv .  (D.14) 

The  discrete  data  is  defined  this  way  until  the  first  underlying 
displacement  dA  s. 

The  slope  of  the  linear  region  intercepts  the  abscissas  at  a  point  dO. 
The  slope  of  the  linear  region  is  the  bundle  modulus  E  (see  Figure  20. 

E  =  (Pui  -  P.„)/(dA1  -  delLJ  (D.15) 

The  slope  can  be  solve  using  £ 

dO  =  delL,  -  P,  „/E  (D.16) 
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LOAD  CURVE 


Figure  2C.  Linear  calculations.  J 

Subsequent  discrete  data  points  can  be  created  by  a  regressive 
process.  At  each  fiber  loss  a  new  slope  emanates  radial  outward  from  dO. 

Between  breaks,  to  define  the  discrete  data,  the  slope  is  defined  by  the  failure  j 

load  Pu,  the  displacement  dA  and  dO 

E  =  [Pu/(dAj  -  dO)]  (&,  -dO)  (D.17) 

where:  dA,  dAj;j=l,2, . n 


6.  Noise 

Noise  was  simulated  simply  by  setting  a  tolerance  level  and  using  the 
random  number  generator  to  produce  a  set  of  random  numbers  between  -0.5 
and  0.5.  This  data  was  added  to  each  discrete  load  £. 


65 


The  input  file  to  the  simulation  has  units  of  kg, mm  and  the 
experimental  output  from  the  Instron  Software  is  in  lbm.in.  A  conversion  of  the 
data  is  required,  this  is  accomplished  through  a  separate  program 

B.  SIMULATION  INTERPRETATION 

The  simulation  was  written  in  fortran  code,  listed  in  Appendix  F.  Multiple 
simulations  were  run  varying  the  parameters  slack  and  noise  to  see  the  effect. 

Simulation  data  was  analyzed  via  the  algorithm  as  represented  in 
Appendix  C.  The  output  is  the  bundle  modulus  E,  the  shape  parameter  of 
distribution  Alpha,  and  the  location  parameter  Beta  (mm/mm). 

The  input  parameters  were  alpha  =  4.0,  beta  =  0.16  kg  and  a  50.0  mm 
Gage  length  (GL).  The  noise  was  varied  from  zero  to  about  0.1%  fluctuation  of 
the  maximum  load.  The  maximum  slack  input  was  3%  of  the  gage  length  and 
a  normal  slack  distribution.  The  seed  used  for  the  random  number  generator 
was  held  constant  so  the  underlying  strength  distribution  was  constant.  The 
seed  for  the  slack  region  and  noise  was  not  held  constant.  The  results  of 
multiple  simulation  test  is  in  Table  IC.  A  value  of  4.134  for  the  shape 
parameter  was  retrieved  consistently  (same  seed),  as  expect  no  deviation  for  the 
same  data.  The  increase  in  slack  drove  the  parameter  down.  The  noise  did  not 
have  any  effect  on  the  parameters.  Figure  3C  show  the  loading  curve  with 
slack  (minus  dO)  and  no  slack. 

No  noticeable  change  in  the  parameter  output  nor  in  the  ECDF  plot 
(Figure  4C). 
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Table  IC.  SIMULATION  OBSERVATIONS  AND 
PROOF  TESTING  RESULTS. 


Maximum  Slack  alpha 
o%  GL 
(B  =  o/2) 

0  4.134 

5  4.010 

10  4.099 

4.095 
4.091 
4.089 
4.088 
4.075 
4.076 


The  numerical  derivative  of  the  slack  region  is  very  sensitive  to  the  any 
noise.  This  can  be  seen  in  the  plot  of  the  slack  ECDF  Figure  5C.  In  order  to 
retrieve  the  slack  distribution  the  recording  system  must  have  a  low  tolerance 
load  cell. 


beta 

(kg) 

0.016 

0.016 

0.016 
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Figure  3C.  Simulation  load  curve,  slack;  no  slack. 
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DISPLACEMENT 


Figure  5C.  Simulation  of  the  slack  ECDF.  Weibull  approximation  of  Normal. 
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APPENDIX  E 


f.  DATA  ANALYSIS  CODE 
A.  ANALYSIS  PROGRAM 

1.  Data  Input 

The  program  requires  a  formatted  data  file  called  EXPER.OUT.  The 
data  is  read  in  and  prompts  the  user  for  several  variables.  These  are  listed 
below. 

-  What  is  the  gage  length? 

-  How  many  fibers  in  the  bundle? 

-  Which  Numerical  differentiation  is  desired? 

3,5,7  pt.  forward  diflf.  or  3,5  pt.  central  diflf. 

-  The  end  of  the  slack  region  and  the  beginning  of  the 
failure  region  in  displacement  (mm)  (this  requires  a  graphical 
estimate). 

-  How  many  points  in  the  input  data  file? 

-  The  acquisition  tolerance  of  the  displacement?  (this  is  for 
any  filtering  of  data  points  in  the  event  the  sampling  rate 
is  too  high  and  repeated  data  is  evidenced. 

-  The  modulus  is  calculated  and  displayed  and  then  ask  if  this 
value  needs  to  be  change.  This  is  in  the  event  the  modulus 
calculation  in  correct  a  manual  entry  can  be  made.  (Be  Careful 
not  to  corrupt  the  true  data). 

2.  Data  Output 


The  data  out  put  is  to  the  screen  and  to  a  file. 


i 


list  of  the  file  outputs. 


The  following  is  a 


-  EXPER.OUT,  the  input  file  is  output  with  an  added  column  of 
strain 

*  SLACK.OUT,  the  slack  region  (load  vs  displacement) 


I 


71 


-  EXMS.OUT,  the  load  curve  minus  slack  region 

-  UBFECDF.OUT,  the  Upper  bound  ECDF 

-  FECDF.OUT,  the  ECDF 

-  SECDF.OUT,  the  slack  region  ECDF 

B.  MAXIMUM  LIKELIHOOD  ESTIMATOR 

The  MLE  program  reads  one  of  the  ECDF  data  files,  users  choice,  the  first 
line  in  the  input  file  is  a  utility  line  defining  the  file  size,  the  number  of 
filaments  and  the  bundle  modulus  E,  the  remaining  lines  are  data. 

1.  Data  Input 

The  input  other  than  the  read  file  is  prompted  from  the  screen.  The 
following  is  a  list  of  requested  information. 

-  How  many  points  for  the  Expected  Ranking  (The  MLE  requires  a 
probability  ranking) 

-  First  approximation  for  alpha 

2.  Data  Output 

The  parameter  estimations  alpha  and  beta  (mm/mm)  are  output  to  the 

screen. 

C.  LINEARIZED  PLOT 

This  program  takes  the  ECDF  and  does  a  Weibull  linearizing  (equation 
C.20). 
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II.  FORTRAN  CODE 
A.  ANALYSIS  PROGRAM 


PROGRAM  ANALY 

*  BUNDLE  FAILURE  ANALYSIS 

* 

* 

*  Outline  of  program: 

*  Read  in  the  data  points 

*  NLA  -  Number  In  Array 
************************************************************** 

*  Find  the 

*  -linear  region 

*  -slope  of  the  linear  region  =  EB  (E  BAR) 

*  -intercept  point  of  the  linear  region  with 

*  the  abscissas  =  DO  (d  zero) 

*  -Maximum  slack  or  minimum  linear  displacement  -  DPI 

*  -Maximum  displacement  of  the  linear  region  -  DP2 
************************************************************** 


*  Input  the  range  of  the  linear  region  by  manually  retrieving 

*  the  defining  points  DPI  and  DP2. 

* 

************************************************************ 

*  OUTPUT:  SLACK. OUT  -  slack  region. 

*  SECDF.OUT  -  slack  ECDF. 

*  EXMS.OUT  -  experimental  data  minus  slack  region. 

*  FECDF.OUT  -  failure  region  ECDF. 

*  UBFECDF.OUT  -  upper  bound  of  the  failure  region  ECDF. 

* 

INTEGER  NIA.NISA 
PARAMETER(NIA=6000NISA=1000) 

REAL  D(NIA),P(NIA),PP(NLA),H,LP(NIA),LD(NIA),FS(NIA),F(NIA) 
+,FUB(NIA),DF(NIA),DFUB(NIA),PF(NIA),PFUB(NIA),TEMP(NIA),GL 
INTEGER  COUNT, NIDP1,NIDP2,NI(NIA),SLCT 
CHARACTERS,  Q1,Q2 

* 

***************************************************************** 

*  Input  analysis  control  parameters. 

* 

***************************************************************** 

PRINT  *,  *What  is  the  gage  length?  (mm)’ 

READ  *,  GL 


PRINT  *,  ’  How  many  filements?’ 

READ  *,  N 

PRINT  *,  ’  Which  method  do  you  require  for  reduction  of 
PRINT  *,  ’  the  slack  region  to  the  ecdf  space?’ 
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PRINT  *,  ’  1)  3  PT  forward  diff.  method.’ 

PRINT  *,  ’  2)  5  PT  forward  diff.  method.’ 

PRINT  ’  3)  7  PT  forward  diff.  method.’ 

PRINT  ’  4)  3  PT  central  diff.  method.’ 

PRINT  ’  5)  5  PT  central  diff.  method.’ 

READ  *,  SLOT 

* 

***************************************************************** 

*  Read  the  range  of  the  linear  displacement. 

* 

***************************************************************** 


PRINT  *,  ’INPUT  DPI  (the  maximum  slack)’ 

READ  *,  DPI 

PRINT  *,  ’INPUT  DP2  (the  maximum  linear  displacement)’ 
READ  *,  DP2 

******************************************************************** 

*  Find  the  how  many  points  are  in  the  data  file. 

* 

******************************************************************** 

PRINT*,  ’How  many  points  in  the  data  file?’ 

READ  *,  COUNT 

******************************************************************** 

*  Find  out  what  the  tolerance  of  the  machine  is. 

*  (displacement) 

******************************************************************** 

PRINT  *,  What  is  the  tolerance  of  the  displacement  output’ 
PRINT  *,  ’data?  This  is  the  limit  between  data  displacement’ 
PRINT  *,  ’Input  in  mm.’ 

READ  *  TOLD 

PRINT  *,  ’How  many  points  in  the  empirical  rank?’ 

READ  *  NER 

144  PRINT  *,  T)o  you  want  the  load  in  Kg  (K)  or  Newtons  (N)?’ 
READ(*,303)  Q2 

IF(Q2  .EQ.  TC  .OR.  Q2  .EQ.  ’N’)  THEN 
PRINT  *,  ’  ’ 

ELSE 

PRINT  *,  Try  again  in  UPPER  CASE  K  OR  N’ 

GOTO  144 
END  IF 


******************************************************************** 

*  Read  in  the  experimental  results. 

* 

******************************************************************** 


OPEN(UNIT=  1 1  ,FILE='EXPER.  OUT,STATUS=’OLD’) 

* 

DO  10  1=1, COUNT 
READU1.110)  NI(I),TEMP(I),P(I) 

10  CONTINUE 
REWINDUl) 
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CONVER  =  1.0 
IF(Q2  .EQ.  *N)  THEN 
CONVER  =  9.8062 
END  IF 

DO  44  1=1, COUNT 
P(I)  =  P(I)*CONVER 

WRITE(ll.lOOl)  NI(I),TEMP(I),P(I),TEMP(iyGL 
44  CONTINUE 

******************************************************************** 

*  Filter  out  any  points  with  difference  in  displacement  interval 

*  less  than  the  machine  tolerance  TOLD. 

* 

******************************************************************** 

J=1 

D(l)  =  TEMP(l) 

DO  11  1=2, COUNT 
DIFF  =  TEMPd)-TEMP(I-l) 

IF(DIFF  .LT.  TOLD)  THEN 
PRINT  *,  ’SKIPPED  POINT  ’,1 
ELSE 
J=J+1 

D(J)  =  TEMP(I) 

P(J)  =  P(I) 

END  IF 

11  CONTINUE 
COUNT  =  J 

* 

******** 

*  FIND  COUNT  (NI)  OF  DPI  AND  DP2 
******** 


1=1 

20  IF(D(I)  .LT.  DPI)  THEN 

1=1+1 
GOTO  20 
ELSE 
NIDP1  =  I 
DPI  =  D(I) 

ENDIF 

30  IFdXI)  .LT.  DP2)  THEN 

1=1+1 
GOTO  30 
ELSE 

NIDP2  =  1-1 
DP2  =  D(I-1) 

ENDIF 

* 

******************************************************************** 
*  OUTPUT  THE  SLACK  REGION. 


******************************************************************** 


0PEN(UNIT=13,FILE=’SLACK0UT,STATUS=TJNKN0WN’) 

DO  500  I  =  1.NIDP1 
WRITE(13,110)  I,D(I),P(I) 

500  CONTINUE 
CLOSE(13) 

* 

******************************************************************** 

*  LEAST  SQUARES  TO  DETERMINE  THE  SLOPE  BETWEEN  DPI  AND 

DP2. 

* 

******************************************************************** 


TOP  =  0.0 
BOT  =  0.0 

DO  40  I=(NIDP1+1),NIDP2 
LP(I)  =  P(I)  -  P(NIDPl) 

LD(I)  =  IXI)  -  DPI 
TOP  =  TOP  +  LP(I)*LD(I) 

BOT  =  BOT  +  LD(I)**2 
40  CONTINUE 

EB  =  TOP/BOT 

********* 
*  Find  the  intercept  point  do 

********* 


I 


303 

* 


DO  =  DPI  -  P(NIDP1)/EB 
PRINT  *,  ’Do  =  ’,D0 
PRINT  *,  ’  ’ 

PRINT  *,  The  Modulus  calculated  from  the  load  curve’ 
PRINT  *,  ’is  ’,EB 
PRINT  *,  ’  ’ 

PRINT  *,  Would  you  like  to  change  this  (Y/N)?’ 
READC^.303)  Q1 

IF(Q1  .EQ.  ’Y’  .OR.  Q1  .EQ.  >’)  THEN 
PRINT  *,  ’Input  the  new  modulus’ 

READ  *,  EB 

ENDIF 

FORMATCA1) 


******************************************************************** 

*  OUTPUT  THE  EXPERIMENTAL  DATA  MINUS  SLACK. 

* 

******************************************************************** 

OPEN(UNIT=13,FILE=’EXMS.OUT’,STATUS=TJNKNOWN’) 

J=1 

WRITE(13,110)  J,0.0,0.0 
DO  550  I  =  NIDPl, COUNT 
J=J+1 

WRITE!  13,1 10)  J,(D(I)-D0),P(I) 

550  CONTINUE 
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******************************************************************** 

*  NUMERICAL  DIFFERENTIATION  OF  THE  SLACK  REGION 

*  (using  a  forward  and  central  difference  methods). 

* 

******************************************************************** 


IF(SLCT  .EQ.  1)  GOTO  103 
IFCSLCT  .EQ.  2)  GOTO  105 
IF(SLCT  .EQ.  3)  GOTO  107 
IF(SLCT  .EQ.  4)  GOTO  130 
IF(SLCT  .EQ.  5)  GOTO  150 

* 

********************************** 

*  3  pt  forward  diff.  method. 

********************************** 


103  DO  3  I=1,NTDP1 

*  if  out  of  the  slack  region  leave  the  loop 

IF  (D(I)  .GT.  DPI)  GOTO  51 
H  =  ABS(D(I+1)-D(I)) 

*  3  point  forward  difference  numerical  differentiation. 

*  Generate  Fs(da)  CDF 

PPG)  =  (-3*P(I)+4*P(I+l)-P(I+2))/(2*H) 

FS(I)  =  PP(I)/EB 
3  CONTINUE 
GOTO  51 

********************************** 

*  5  pt  forward  difference  method. 
********************************** 

105  DO  5  I=1,NIDP1 

*  if  out  of  the  slack  region  leave  the  loop 

IF  (D(I)  .GT.  DPI)  GOTO  51 
H  =  ABS(D(I+1)-D(I)) 

*  5  point  forward  difference  numerical  differentiation. 

*  Generate  Fs(da)  CDF 

PPG)  =  (-25*P(I)+48*P(I+l)-36*PG+2)+16*PG+3)- 
C3*P(I+4))/(  12*H) 

FS(I)  =  PP(I)/EB 

5  CONTINUE 
********************************** 

*  7  pt  forward  difference  method 
********************************** 

107  DO  7  I=1,NIDP1 

*  if  out  of  the  slack  region  leave  the  loop 

IF  (D(I)  .GT.  DPI)  GOTO  51 
H  =  ABSGXI+D-DG)) 

*  7  point  forward  difference  numerical  differentiation. 

*  Generate  Fs(da)  CDF 

PPG)  =  (-147*P(I)+360*P(I+l)-450*PG+2)+ 
+400*P(l+3)-225*PG+4)+72*PG+5>-10*PG+6))/(60*H) 
FS(I)  =  PP(I)/EB 
7  CONTINUE 
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********************************** 

*  3  pt  central  difference  method 
********************************** 


* 

130  1=1 

H  =  ABS(DG+1)*DG)) 

PP(I)  =  (-3*P(I  )+4*P(I+ 1  )*P(I+3 )  )/(2  *H) 
FS(I)  =  PPdVEB 
DO  33  I=2,NIDP1 

*  if  out  of  the  slack  region  leave  the  loop 

IF  (D(I)  .GT.  DPI)  GOTO  51 
H  =  ABS(D(I+1)-D(I)) 

PP(I)  =  (-P(M>+P(I+l)y(2*H) 

*  Generate  Fs(da)  CDF 

FS(I)  =  PP(IVEB 

33  CONTINUE 
********************************** 

*  5  pt  central  difference  method 
********************************** 


* 

150  1=1 

*  Use  a  3  pt  forward  difference  method  for  first  couple  pts. 

H  =  ABSdXI+l)-D{I)) 

PPd)  =  (-3*Pd)+4*PO+l)-Pd+3))/(2*H) 

FSd)  =  PP(iyEB 
1=2 

H  =  ABS(Dd+l)-Dd)) 

PPd)  =  (-3*P(D+4*P(I+l)-P(I+3))/(2*H) 

FS(I)  =  PPdVEB 
DO  55  I=3,NIDP1 

*  if  out  of  the  slack  region  leave  the  loop 

IF  (Dd)  -GT.  DPI)  GOTO  51 
H  =  ABSdXI+D-Dd)) 

PPd)  =  (P(I-2)-8*P(I-l  )+8*P(I+l  )-P(I+2))/(  12*H) 

*  Generate  Fs(da)  CDF 

FSd)  =  PPdVEB 

55  CONTINUE 
********************************** 

51  CONTINUE 

******************************************************************** 

*  FAILURE  REGION 

* 

******************************************************************** 

J=0 

DO  60  I=NTDP2, COUNT 
J=J+1 

F(J)  =  1  -  PdV(EB*(Dd)-D0)) 

DF(J)  =  DdVDO 
60  CONTINUE 
NUMF  =  J 

******************************************************************** 
*  EXTRACT  THE  UPPER  BOUND. 
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* 

******************************************************************** 
DEL  =  1.0/NER 
K=NUMF 
I=NER 

FUBG+l)  =  F(K) 

DFUBG+l)  =  DF(K) 

K=K-1 

1=1-1 

FL  =  1.0*I/NER  -  DEL/2 
FU  =  1.0*I/NER  +  DEL/2 
PRINT  *,TL  F(K)  FU’ 

PRINT  *  FL,’  <  '?(K);  <  \FU 

90  IF((K  .LT.  1)  .OR.  (I  .LT.  0))  GOTO  92 

91  IF((F(K)  .LT.  FU)  .AND.  (F(K)  .GE.  FL))  THEN 

FUBG+l)  =  F(K) 

DFUBG+l)  =  DF(K) 

1=1-1 

FL  =  1.0*I/NER  -  DEL/2 
FU  =  1.0G/NER  +  DEL/2 
K=K-1 

IF(K  .LT.  1)  GOTO  92 
GOTO  91 

ELSEIF(F(K)  .LT.  FL)  THEN 

******** 

*  check  to  see  how  many  it  skips 

******** 

J1  =  1 

93  IF(F(K)  .LT.  (1.0*G-J1)/NER  -  DEL/2))  THEN 

J1=J1+1 
GOTO  93 
ENDIF 

DO  94  I1=I,(I-J1+1),-1 
IFG1+1  .LT.  1)  GOTO  92 
FUBG1+1)  =  1.0*I1/NER 
DFUBG1+1)  =  DF(K) 

94  CONTINUE 

1=11 

FUBG+l  )=F(K) 

DFUB(I+1  )=DF(K) 

1=1-1 

FL  =  1.0*I/NER  -  DEL/2 
FU  =  1.0*I/NER  +  DEL/2 
K=K-1 
GOTO  90 
ELSE 
K=K-1 

IF(K  .LT.  1)  GOTO  92 
GOTO  90 
ENDIF 

92  CONTINUE 

PRINT  *  ’OUT  OF  LOOP  WITH  K  AND  I  ’ 
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PRINT  *,  'K=  \K,'  1=  M 
NUMEF  =  1+2 

************************************************************** 

*  WRITE  THE  DATA 

* 

************************************************************** 

OPEN  (UNIT=ll,FILE=TJBFECDF.OUT’,STATUS=’UNKNOWN’) 
MSUM  =  NER+1 -NUMEF 
WRITE(  11,1010)  MSUM.N.EB 

* 

DO  250  I=NUMEF,NER 
WRITE(11,110)  I,DFUB( I ),FUB( I ) 

250  CONTINUE 
CLOSE(ll) 

OPEN(UNIT=ll,FILE=’FECDF.OUT’,STATUS='UNKNOWN’) 
WRITE(11,1010)  NUMF,N,EB 
DO  70  I-1.NUMF 
WRITE(ll.llO)  I,DF(I),F(I) 

70  CONTINUE 
CLOSE(ll) 

OPEN(UNIT=12,FILE=’SECDFUUT,STATUS=,UNKNOWN’) 

DO  80  I=1,NIDP1 
WRITE(12,110)  I,Dd),FS(I) 

80  CONTINUE 

*  *  4c  4c  *  *  4c  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  ♦  *  *  *  *  * *  *  * *  *  * 

*  FORMATS 

* 

******************************************************************** 

100  F0RMAT(1X,I5) 

110  F0RMAT(1X,I5,2X,F8.4,2X,F8.4) 

1001  FORMAT(1X,I5,2X,F8.4,2X,F8.4,2X,E10.4) 

1000  FORMAT(  1X,I5,2X,E16. 10,2X,E  16. 10) 

1010  FORMAT(  1X,I5,2X,I5,2X,F8.4) 
************************************************************ 

END 

B.  MAXIMUM  LIKELIHOOD  ESTIMATOR 

PROGRAM  MLE 
INTEGER  NIA 
PARAMETER  (NIA  =  6000) 

REAL  FE(NIA),DFE(NIA),FUB(NIA),DFUB(NIA),GL 
CHARACTER*  1  Ql 

* 

PRINT  *,  ’Do  you  want  the  upper  bound  or  complete  data’ 
PRINT  *,  ’set,  suggest  the  compete  data  set  when  the’ 

PRINT  *,  ’DATA  is  less  the  the  number  of  filements’ 

110  PRINT  *,  TYPE  "U"  for  upper  and  "C”  otherwise’ 

303  FORMATXAl) 

READ  (*,303)  Ql 
IF(Q1  .EQ.  ’U)  THEN 

OPEN(UNIT=  1 1  ,FILE=UB  FECDF.  OUT,STATUS=’OLD’) 
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k 


I 


I 


I 


■ 


I 

» 


ELSEIF  (Ql  .EQ.  ’C’)  THEN 

OPEN(UNIT=  1 1  tFILE=’FECDF.OUT,STATUS=’OLD’) 

ELSE 

PRINT  *  TRY  TYPING  CAPITAL  LETTERS’ 

GOTO  110 
ENDIF 

* 

READdl.1010)  M,N,EB 
1010  F0RMAT(1X,I5,2X,I5,2X,F8.4) 

DO  10  I=1^M 

REAIX1 1,1020)  J,DFUB(I),FUB(I) 

10  CONTINUE 

1020  F0RMAT(1X,I5,2X,F8.4,2X,F8.4) 

************************************************************** 

*  CALCULATE  ALPHA  AND  BETA 

PRINT  *,  'What  is  the  GAGE  LENGTH?’ 

READ  *,  GL 

* 

*  REDUCE  THE  DATA  TO  ABOUT  50  POINTS 
888  PRINT  *,  How  many  points  for  the  MLE  ?’ 

READ  *,  MLE 
PRINT  *,  MLE 
DEL  =  1.0/(MLE+1.0) 

J  =  1 
I  =1 

*  find  the  value  of  F  above  and  below  the  expected  rank  and  average 
400  FE(I)  =  1.0*I/(MLE+1.0) 

IF(I  .GT.  MLE)  GOTO  450 

410  IF((FUB(J)  .LE.  FE(I))  .AND.  (FUB(J+1)  .GT.  FE(I)))  THEN 

DFE(I)=  -(DFUB(  J+ 1  )-DFUB(  J))*(FUB(  J+ 1  )-FE(I))/(FUB(  J+ 1  )-FUB(  J)) 
+  +  DFUB(J+1) 

1=1+1 

GOTO  400 

ELSEIF((J+1)  .LE.  (M))  THEN 

J=J+1 

GOTO  410 

ENDIF 

450  PRINT  *,  I, J,FUB( J),FUB( J + 1),FE(I ) 

IF(FUB(J+1)  .LE.  0.0)  THEN 
PRINT  *,  ’MLE  to  large  try  again’ 

GOTO  888 
ENDIF 

************************************************************** 

*  estimate  the  initial  alpha 

PRINT  *,  ’Estimate  the  intial  alpha’ 

READ  *  ALFA0 

************************************************************** 

*  calculate 

ALFA  =  ALFA0 
KK=0 

320  DO  310  IT  =  1,2 
SUM1=0.0 
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SUM2=0.0 

SUM3=0.0 

*  Calculate  the  variables  of  alpha 

DO  300  1=1, MLE 
IF(DFE(I)  .GT.  0.0)  THEN 

SUM1  =  SUM1  +  ((DFE(I))**ALFA  *  LOG{DFE(I))) 

SUM2  =  SUM2  +  ((DFE(I))**ALFA) 

SUM3  =  SUM3  +  LOCKDFE(I)) 

END  IF 

300  CONTINUE 

EF(IT  .EQ.  1)  THEN 

ALFA1=1.0/(SUM1/SUM2  -  (SUM3/MLE)) 

ALFA  =  ALFA1 
ELSEIF(IT  .EQ.  2)  THEN 
ALFA2= 1 .0/(SUM  1/SUM2  -  (SUM3/MLE)) 

ENDIF 

310  CONTINUE 

ALFA0P  =  ALFA0  -((ALFA0  -  ALFA  1  )*  *2  >'(  ALFAO  -  2*ALFA1 
+  +  ALFA2) 

TEST  =  ABS(ALFAOP-ALFAO) 

IF(TEST  .GT.  .0001)  THEN 
IF(KK  .GT.  1000)  GOTO  341 
KK  =  KK  +  1 
ALFAO  =  ALFAOP 
PRINT  *  ’ALFAO  =’,ALFA0 
GOTO  320 
ELSE 

ALPHA  =  ALFAOP 

ENDIF 

SUM  =  0.0 

DO  340  I  =  l.MLE 

SUM  =  SUM  +  (DFE(I))**ALPHA 

340  CONTINUE 

341  BETA  =  (SUM/(MLE))**(  1/ALPHA) 

PRINT  *,  ’ALPHA  =  ’.ALPHA,’  BETA  =  ’.BETA/GL,’  mm/mm’ 
END 

C.  WEIBULL  UNEAR  CDF 

PROGRAM  LIN 
INTEGER  NLA, AREA 
PARAMETER  (NIA  =  6000) 

REAL  FE(NIA),DFE(NIA),FUB(NIA),DFUB(NIA),FS(NIA),D(NIA) 

* 

0PEN(UNIT=11,FILE=TJBFECDF.0UT,STATUS=’0LD’) 
READ(11,1010)  M,N,EB 
1010  F0RMAT(1X,I5,2X,I5,2X,F8.4) 

DO  10  1=1, M 

READ(1 1,1020)  J,DFUB(I),FUB(I) 

10  CONTINUE 

1020  F0RMAT(1X,I5,2X,F8.4,2X,F8.4) 

************************************************************** 
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*  REDUCE  THE  DATA  TO  ABOUT  50  POINTS,  USING  THE  EXPECTED 
RANK 

PRINT  *,  HOW  MANY  POINTS  TO  FILTER  USING  THE  EXPECTED 

RANK 

READ  *  NUM 
PRINT  *  NUM 
DEL  =  1.0/(NUM+1.0) 

J  =  1 
I  =1 

*  find  the  value  of  F  above  and  below  the  expected  rank  and  average 
400  FE(I)  =  1.0*I/(NUM+1.0) 

IF(I  .GT.  NUM)  GOTO  450 

410  IF((FUB(J)  .LE.  FE(I))  .AND.  (FUB(J+1)  .GT.  FE(I)))  THEN 

DFE(I)=  -(DFUB(J+1)-DFUB(J))*(FUB(J+1)-FE(I))/(FUB(J+1)-FUB(J)) 

+  +  DFUB(J+1) 

1=1+1 

GOTO  400 

ELSEIF((J+1)  .LT.  (M+l))  THEN 

J=J+1 

GOTO  410 

ENDIF 

450  PRINT  *  I,J 

******************************************************************** 

*  CALCULATE  ALPHA  AND  BETA 

* 

*  FIND  F*=LOG<-LOG(l-F))  VS  LOG{DELT A/BETA) 

DO  20  1=1, M 

IF(FUB(I)  .LE.  0.0)  FUB(I)  =  0.0001 
FS(I)  =  LOG(-LOGU-FUBG))) 

IF(DFUB(I)  .LE.  0.0)  DFUB(I)  =  0.0001 
D(I)  =  LOG<DFUB(I)) 

20  CONTINUE 

DO  70  I=1,M 

510  IF((FS(I)  .LE.  0.63212)  .AND.  (FSd+1)  -GT.  0.63212))  THEN 
BETA  =  -(D(I+l)-Dd))*(FS(I+l)-0.63212y(FSd+l)-FSd)) 

+  +  Dd+1) 

GOTO  71 
ENDIF 

70  CONTINUE 

71  CONTINUE 
BETA  =  EXP(BETA) 

PRINT  *,  ’  BETA  =  ’.BETA 


DO  80  1=1, M 

IF(DFUB(I)  .LE.  0.0)  DFUBd)  =  0.0001 
D(I)  =  LOGQFUBdVBETA) 

80  CONTINUE 

BETA  =  EB*BETA/N 

************************************************************** ****** 

OPEN(UNIT=ll,FILE=’FSTAR.OUT,STATUS=’UNKNOWN’) 

DO  60  1=1, M 
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WHITE(11,1020)  I,D(I),FS(I) 

CONTINUE 

END 
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APPENDIX  F 


SIMULATION  FORTRAN  CODE 
PROGRAM  SIMUL 


SIMUL  IS  A  STATISTICAL  SIMULATION  ROUTINE  USED 
FOR  RELIABILITY 

ANALYSIS  OF  COMPOSITE  FILEMENTS. 


UTILITIES  FILES  REQUIRED  ARE: 

SEED  -  file  containing  the  integer  seed  for  random 
number  generation 

INPUT  -  file  containing  the  input  parameters  * 


OUTPUT  FILES  INCLUDE: 

EXPER.OUT  -  the  experimental  simulation  of  discrete 
data. 

SEED  -  outputs  a  new  seed. 


INPUT  FILE  FOR  BUNDLE2.FOR  PROGRAM  -.line  1 
XXXXXXXXXXXXXXXXXXX12345678  :line  2 

:line  3 
.line  4 
:line  5 
dine  6 
dine  7 
.‘line  8 

:line  9 
dine  10 

dine  11  %  OF  MAX 
DISPLACEMENT 
:line  12  pts/s 
dine  13  mm/min 
-.line  14  tolerance 
:line  15  tolerance  of  displ. 
output 


INPUT  PARAMETERS  ARE: 

1)  THE  LENGTH  OF  THE  TEST  BUNDLE  (MGL). 

2)  THE  NUMBER  OF  FILEMENTS  PER  BUNDLE  (N). 

3)  THE  NUMBER  OF  SEGMENTS  PER  FILEMENT  (M). 

4)  THE  SEGMENT  LENGTH  IS  DEFINED  BY  MGL=M*LS 

THIS  PROGRAM  CREATES: 


E  (kg) 

ROWS  =  M 

COLS  =  N 

MGL  (mm) 

ALPHA 

BETA 

MODEL 

SKEW/ALPHA 


.90 


100 

100 

250.0 

4.0 

16.0 

NORMAL 

3.5 


MAX  SLACK  (%  )  10.0 


SAMPLE  RATE  20.0 

CROSS-HEAD  SPEED  1.8 
NOISE  +/-  kg  .001 
TOLD  (mm)  .00254 


1)  THE  RANDOM  NUMBERS 
‘  2)  CONVERTS  TO  WIEBULL  BUNDLE 

3)  FIND  THE  WEAKEST  IN  THE  ENTIRE  BUNDLE 
FOR  THE  LOWER  BOUND. 

4)  FIND  THE  WEAKEST  IN  EACH  FILEMENT 

5)  ORDER  THE  WEAKEST  TO  STRONGEST  FILEMENTS 

6)  CHANGE  THE  STRENGTHS  TO  BREAKING  LOADS  AND 
DISPLACEMENTS 

************************************************************** 
Begin  Program: 

:************************************************************** 


INTEGER  ROWS, COLS, MAX, LIM, LIMl  ,N PTS, COUNT, LIM2 
PARAMETER(LIM=6000,LIM1=1000,LIM2=1001) 

REAL  BUN(LIM1),RAND1(LIM1),RAND2(LIM1),DISPL(LIM1),  -> 

+SLACK(LIM1),P(LIM1),ALPHA,BETA,MGL,BRKLDM,  * 

+DISPLM,MAXX,MAXY,MEAN,PU(LIM1 ), 
+PL(LIM1),SUM,MAXSLK,TEMP(LIM2),PS(LIM2),D(LIM), 
+DISPLA(LIMl),DEL,PMS,LOAD(LIM)NSE 
DIMENSION  N0(LIM1),NB(LIM1),NDB(LIM1) 

DOUBLE  PRECISION  SEED 

CHARACTER*  1  MODEL*7,SONS,Q,Q2,Q3,Q4  ■ 

COMMON  MAX 
COMMON  /SMIVM,N 

* 

******************************************************************** 


*  FILE  11  =  PARAM  DATA  FOR  INPUT  ■ 

*  FILE  12  =  SEED  DATA  FOR  SEEDING 

* 

******************************************************************** 


OPEN(UNIT=  1 1  ,FILE=TNPUT\STATUS=’OLD’) 

OPEN(UNIT=12,FILE=’SEED',STATUS='OLD’)  I 

* 

************************************************************** 

*  READ  IN  THE  INPUT  DATA 

* 

************************************************************** 


READ(11,1530) 

READ(11,1500)  E 
PRINT  *,  ’MODULUS^,  E 
READd  1,1520)  M 
PRINT  *,TM  =  \M 
ROWS=M 
READdl.1520)  N 
COLS=N 

PRINT  *,  N  =  \N 

READ(11,1510)  MGL 

PRINT*,  ’Mean  Gage  Length  =  ’,MGL 


i 
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MAX=M*N 

REAEK  11,1510)  ALPHA 
PRINT  *,  'ALPHA  =  ’ALPHA 
READ(11,1510)  BETA 
PRINT  *  'BETA  =’,BETA 

* 

****** 

*  Determine  which  model  to  nee  for  slack. 
****** 


READ(11,1515)  MODEL 

PRINT  *  MODELs’MODEL,’  APPROXIMATION’ 

IF  (MODEL  .EQ.  ’NORMAL’)  THEN 
READd  1,1500)  STDDEV 

PRINT  *,  'ALPHA  =  STDDEV 
READ(11,1510)  MEAN 

PRINT  *,  'Maximum  slack  as  percent  of  expected’ 

PRINT  *,  ’maximum  displacement  ( 3 %  gage  length). =’, ME  AN 
MEAN  =  (MEAN/100)*(0.03)*MGL/2 
PRINT  *  'Means’, ME  AN, ’mm’ 

PRINT  *,  Maximum  slack  =’,MEAN*2,’mm’ 

ELSE 

READdl.1500)  RANGE  1 
READ(11,1500)  RANGE2 
RANGE  1=(RANGE  1/100)*(.03)*MGL 
RANGE2=(RANGE2/100)*(.03)*MGL 
PRINT  *,  TtANGE  1  =’, RANGE  1,’  mm’ 

PRINT  *,  -RANGE  2  =’,RANGE2,’  mm’ 

ENDIF 

* 

****** 

*  Read  in  the  sampling  rate  and  cross-head  speed. 
****** 

*************************************** ***** *********************** 

*  Ensure  the  displacement  interval  governed  by  the  cross-head  and 

*  sampling  rate  is  greater  than  the  tolerance  level  of  the  output 

*  displacement  data. 

* 

******************************************************************* 


READC11.1510)  SR 
PRINT  *,  'Sampling  rate  =’,SR 
READU1.1510)  XSP 
PRINT  *,  ’Cross-head  speed  =’,XSP 

****** 
*  How  much  noise. 

****** 


READC11.1510)  NSE 
PRINT  *,  -Range  of  noise  =’,NSE 
READdl.1510)  TOLD 


PRINT  *,  Tolerance  of  output  =’,TOLD 

*  Calculate  the  number  of  total  points  recorded  in  the  * 

*  experiment.  * 

********* 

* 

444  RR  =  MGL*(0.03) 

PTS  =  RR*SR/(XSP/60) 

DELTAD  =  RR/PTS 
IF(TOLD  .GT.  DELTAD)  THEN 

PRINT  *,  The  displacement  interval  is  less  than  the  output* 
PRINT  *,  ’tolerance,  please  input  a  new  cross-head  and/or’ 
PRINT  *,  ’sampling  rate.’ 

PRINT  *,  ’Current  Cross-head  =  ’,XSP 
PRINT  •  Input  new  XSP’ 

READ  *,  XSP 

PRINT  *,  ’Current  sampling  rate  =  ’,SR 
PRINT  *,  ’Input  new  sampling  rate’ 

READ  *,  SR 
GOTO  444 
ENDIF 

* 

****** 

*  Determine  the  type  of  simulation. 
****** 

* 

PRINT  *,  ’SLACK  OR  NO  SLACK  ?  (Y/N)’ 

READ  (*,1600)  SONS 

PRINT  *,  ’NOISE  OR  NO  NOISE  ?  (Y/N)* 

READ  (*,1600)  Q2 

PRINT  *,  ’Change  the  SEED  ?’ 

READ  (*,1600)  Q 

PRINT  *,  ’COMPLAINCE  (Y/N)?’ 

READ  (*,1600)  Q3 

IF((Q3  .EQ.  Y)  .OR.  (Q3  .EQ.  ’y’))  THEN 

PRINT  *,  'WHAT  IS  THE  COMPLIANCE  NOISE  LEVEL?’ 

READ  *,  CTOL 

ENDIF 

PRINT  *,  T)o  you  want  English  units  (lbs, inches)  on  output? 
+(Y/N )’ 

PRINT  *,  ’If  N  the  output  is  in  SI  units  (kg, mm)’ 

READ  (*,1600)  Q4 

* 

CLOSE(ll) 
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******************************************************************** 


*  BEGIN  THE  ALGORITHM 

******************************************************************** 

** 

********* 

*  Read  the  seed,  (seed  for  the  uniform  random  number  generator). 

*  and  generate  the  random  numbers 
********* 

* 

READU2.1000)  SEED 
REWIND*  12) 

DO  10  1=1, N 

CALL  RAND(SEED,M,BUN) 

CALL  SMALL* ROWS, BUN, RAND  1(D) 

10  CONTINUE 


********* 

*  Check  if  slack  is  requested  in  simulation. 
********* 


IF  (SONS  .EQ.  Y)  THEN 
CALL  RAND(SEED,N,RAND2) 

ENDIF 

* 

********* 
*  initialize  the  order  arrays. 

********* 


DO  25  I=1,N 
SLACK(I)=0.0 
NO(I)=I 
NB(I)=I 
NDB(I)=I 

25  CONTINUE 
* 

********* 

*  Check  if  seed  is  to  be  changed,  if  so  write  the  new  seed. 
********* 


IF  *(Q  .EQ.  V)  .OR.  (Q  .EQ.  V))  THEN 
WRITE*  12,1000)  SEED 
REWIND  (12) 

CLOSE*  12) 

ENDIF 


******************************************************************** 
*  Generate  the  load  and  slack  distributions. 


******************************************************************** 


** 
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CALL  WEIBUIX  ALPHA, BET  A,N,RAND1,P) 

IF  (SONS  .EQ.  V  .OR.  SONS  .EQ.  y>  THEN 
IF  (MODEL  .EQ.  ’NORMAL’)  THEN 

CALL  NORMAL(STDDEV,MEAN,N,RAND2, SLACK, MAXSLK) 
ELSE 

CALL  UNIFM(RANGE1,RANGE2,N,RAND2, SLACK, MAXSLK) 
ENDIF 
SLACK(1)=0.0 
ENDIF 

* 

******************************************************************** 
*  Find  the  displacement  for  each  failure  load  and  slack. 

******************************************************************** 


DO  20  1=1, N 

DISPLG)  =  (P(I)/E)*(SLACK(I)  +  MGL) 
20  CONTINUE 
K=N 

****** 

*  Order  the  "DISPL"  array,  and  carry  "No" 
****** 


CALL  SHELUN,DISPL,NO) 

* 

****** 

*  reorder  the  remaining  arrays,  the  others  are 

*  keyed  to  "DISPL” 

****** 


CALL  SWITCH(N,NO,RANDl) 

CALL  SWITCH(N,NO,RAND2) 

CALL  SWITCH(N ,NO, SLACK) 

CALL  SWITCH(N,NO,P) 

* 

******************************************************************** 

*  Calculate  the  apparent  displacement. 
******************************************************************** 


DO  40  1=1  JV 

DISPLAG)  =  DISPLG)  +  SLACK(I) 

40  CONTINUE 

* 

******* 

*  Reorder  the  arrays  keyed  to  "DISPLA" 
******* 


CALL  SHELUN, DISPLAY) 
CALL  SWITCH(N,NB,RAND1) 
CALL  SWITCH(N,NB,RAND2) 
CALL  SWITCH(N,NB,SLACK) 
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CALL  SWITCH(N ,NB  JDISPL) 

CALL  SWITCH(N,NB,P) 

CALL  SWCH2(N,NB,N0) 

* 

******* 

*  Define  the  ordered  slack  as  "TEMP" 

******* 

* 

DO  45  I=1,N 
TEMP(I)=SLACK(I) 

45  CONTINUE 

IF((SONS  .EQ.  V)  .OR.  (SONS  .EQ.  V))  THEN 
CALL  SHELL(N,TEMP,NDB) 

END  IF 

* 

******************************************************************** 
*  * 

*  Calculate  the  maximum  slack  and  intercept  point.  * 

*  Calculate  the  upper  and  lower  Pi  bundle  load  at  each  displa.  * 

*  * 


******************************************************************** 

*  Find  the  loads  at  the  slack  points. 
******* 


PS(1)=0.0 
DO  70  J=2,N 
CONST=0.0 
DO  75  1=1, (J-l) 

CONST=CONST+E*(TEMP(J)-TEMP(I))/(MGL+TEMP(I)) 

75  CONTINUE 
PS(J)=CONST 
70  CONTINUE 

******************************************************************** 

*  ADD  COMPLIANCE  TO  SLACK 

*  * 

*  DC1  =  mm  * 

******************************************************************** 

IF((Q3  .EQ.  V)  .OR.  (Q3  .EQ.  V))  THEN 

DC1  =  .0356 

Al=  0.0512 

A2=  -0.0217 

POWER  =  0.1975 

A3  =  0.0064 

A4  =  0.0242 

DO  120  I  =  1.N 

IF((TEMP(I)  .LT.  DC1)  .AND.  (PS(I)  .GE.  CTOL))  THEN 
TEMP(I)  =  TEMP(I)  +  Al*(PS(I)**POWER)  +  A2 
ELSEIFXPSd)  .GE.  .01)  THEN 
TEMP(I)  =  TEMP(I)  +  A3*PS(I)  +  A4 

FNTYTF 

120  CONTINUE 
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ENDIF 


♦*****♦****♦**♦****♦*♦**♦**♦**♦*♦♦*♦**♦  ********************  ********* 
*  FAILURE  POINTS. 

******************************************************************** 


DO  80  J=1N 
CONST=0.0 
DO  85  I=JN 

CONST=CONST+E*(DISPLA(J)-SLACK(I))/(MGL+SLACK(I)) 

85  CONTINUE 
PU(J)=CONST 
80  CONTINUE 

******************************************************************** 

*  ADD  COMPLIANCE  TO  FAILURE  REGION 

*  * 

*  DC1  =  mm  * 

******************************************************************** 

IF((Q3  .EQ.  Y)  .OR.  (Q3  .EQ.  y»  THEN 
DO  121  I  =  1,N 

IF((DISPLA(I)  .LT.  DC1)  .AND.  (PU(I)  .GE.  CTOL))  THEN 
DISPLAG)  =  DISPLAG)  +  Al*(PU(I)**POWER)  +  A2 
ELSEIF(PU(I)  .GE.  .01)  THEN 
DISPLAG)  =  DISPLAG)  +  A3*PUG)  +  A4 
ENDIF 

121  CONTINUE 
ENDIF 

* 

******************************************************************* 

*  ADD  DISCRETE  DATA 

*  TO  THE 

*  SLACK  REGION. 

* 

******************************************************************* 


J=1 

ADD  =  0.0 
CONST=0.0 
EKJ)=DELTAD 
DO  90  1=1, (N-l) 

IF  (IXJ)  .LT.  TEMPG+1))  THEN 
SLOPE  =  (PSG+D-PSGMTEMPG+U-TEMPG)) 
91  LOAD(J)=(D(J)-TEMPG))*SLOPE  +  ADD 
CONST=D(J) 

J  =  J  +  1 

D(J)=CONST+DELTAD 
IF  (D(J)  .LT.  TEMPG+1))  GOTO  91 
ENDIF 

ADD  =  ADD  +  PSG+U-PSG) 

90  CONTINUE 
JBIGL  =  J-l 
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******************************************************************* 

*  ADD  DISCRETE  DATA 

*  TO  THE 

*  LINEAR  REGION 

* 

******************************************************************* 

* 

CONST=0.0 

SLOPE=(PU(l)-PS(N)y(DISPLA(l)-TEMP(N)) 

100  LOAD(J)=SLOPE*(D(J)-TEMP(N))  +  ADD 
CONST=  D(J) 

J=J+1 

D(J)=CONST  +  DELTAD 

IF  (D(J)  .LT.  DISPLA(l))  GOTO  100 

JENDL  =  J-l 

* 

****** 

*  Find  the  intercept  point. 
****** 

* 

DINT=  DISPLAd)  -  PU(l)/SLOPE 

* 

******************************************************************* 

*  ADD  DISCRETE  DATA 

*  TO  THE 

*  FAILURE  REGION. 

* 

******************************************************************* 
DO  110  1=2, N 

IF  (EKJ)  .LT.  DISPLAd))  THEN 
SLOPE  =  PU(I)/(DISPLA(I)-DINT) 

115  LOAD(J)=SLOPE*(D(J)-DINT) 

CONST  =  D{J) 

J=J+1 

D(J)  =  CONST  +  DELTAD 
IF  (D(J)  .LE.  DISPLAd))  GOTO  115 
ENDIF 

110  CONTINUE 
LOAD(J)=0.0 

700  PRINT  *,  ’CONTINOUS  DATA  COMPLETE  * 

COUNT  =  J 
PRINT  *,  'COUNTS, J 

* 

******* 

*  Find  the  largest  value  for  plotting 
******* 

* 

CALL  LARGE(K,PU3RKLDM) 

* 

****** 

*  Where  brkldm  and  sldm  are  the  maximum  for  that  vector. 
****** 
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* 

DISPLM=DISPLA(N) 

* 

****************************************************************** 

*  CREATE  THE  NOISE 

* 

****************************************************************** 


* 

IF  «Q2  .EQ.  Y)  OR.  (Q2  .EQ.  y»  THEN 

CALL  NOISE(SEED, COUNT, LOAD,NSE)  I 

ENDIF 

* 

******************************************************************** 

*  PRINT  THE  DATA  * 

*  * 

******************************************************************** 

OPEN(UNIT=13,FILE=’EXPER.OUT,STATUS=TJNKNOWN’) 

IF((Q4  .EQ.  Y)  .OR.  (Q4  .EQ.  y))  THEN 
CONST1=0.03937 
CONST2=2.2046 
ELSE 

CONST1  =  1.0 
CONST2  =  1.0 
ENDIF 

* 

DO  105  1=1, COUNT 

WRITE(13,111)  I,(D(I)*CONSTl),(LOAD(I)*CONST2) 

105  CONTINUE 

WRrrE(13,116)  ALPHA, BETA, N 
WRITE(13,111)  COUNT,DISPLM,BRKLDM 
CLOSE(13) 

******************************************************************** 

*  FORMATS  * 

*  * 

******************************************************************** 

111  F0RMAT(1X,I5,2X,F8.4,2X,F8.4) 

116  FORMAT  ( 1X,’PARAMETERSALPHA=’,F4 . 1 ,2X,’BETA=\ 

+F4.U1X 

+/N=\I4,/,  1X,’N,,5X,’DELTA’)  1  OX, ’LOAD’) 

1000  F0RMAT(1X,F15.1) 

1500  FORMAT(20X,E16.10) 

1510  FORMA1X20X,F8.4) 

1515  FORMAT(20X,A6) 

1520  FORMAT(20X,I5 ) 

1530  FORMATUX/) 

1600  FORMATXAl) 

* 

STOP 
END 

* _ 

SUBROUTINE  NOISE(SEED,COUNT,LOAD,NSE) 


I 
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INTEGER  COUNT 
DOUBLE  PRECISION  SEED 
REAL  LOAD(COUNT)NSE,RNDM 
DO  10  1=1, COUNT 
CALL  RANDN(SEED,RNDM) 

LOAD(I)  =  LOAD(I)  +  (RNDM*NSE  -  NSE/2) 
10  CONTINUE 
END 


SUBROUTINE  SWITCH(N,ORDER,ARRAY) 


INTEGER  N,ORDER(N),LIM  1 
PARAMETEIK  LIM=5000) 
REAL  ARRAY(N),TEMP(LIM) 
DO  1  I=1,N 
TEMP(D=ARRAY(I) 

1  CONTINUE 
* 

DO  2  I=1,N 
K=ORDERO) 
ARRAY(I)=TEMP(K) 

2  CONTINUE 
END 


SUBROUTINE  SWCH2(N,ORDER,ARRAY) 


INTEGER  N.LIMl 
PARAMETEIK  LIM1 =5000) 

INTEGER  ARRAY(N),TEMP(LIMl),ORDER(N) 
DO  1  I=1,N 
TEMP(I)=ARRAYd) 

1  CONTINUE 
* 

DO  2  I=1,N 

ARRAY(I)=TEMP(ORDER(I)) 

2  CONTINUE 
END 
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in  * 


SUBROUTINE  LARGE(MAX,  VECTOR, MAXVAL) 


INTEGER  MAX 

REAL  VECTOR(MAX), MAXVAL 
MAXVAL=VECTOR(  1 ) 

DO  10  J=2,MAX 

IF  (VECTOR(J)  .GT.  MAXVAL)  THEN 
MAXVAL  =  VECTOR(J) 

ENDIF 

10  CONTINUE 
END 

SUBROUTINE  SMALUTOTAL, VECTOR, LOW) 


*  THIS  ROUTINE  TAKES  A  MATRIX  IN  VECTOR  FORM  AND 

*  FINDS  THE  SMALLEST 

*  VALUE  IN  EACH  ROW  AND  PUTS  IT  INTO  A  VECTOR 

INTEGER  TOTAL, ROWS, COLS 
REAL  VECTOR(TOTAL),LOW 
COMMON  /SMIV  ROWS, COLS 

*  LOOK  AT  ONLY  ROWS 

LOW= VECTOR!  1 ) 

* 

DO  10  1=1, ROWS 

IF  (VECTORG)  .LT.  LOW)  LOW=VECTOR(I) 

* 

10  CONTINUE 
END 


SUBROUTINE  SHELIXMAXNUM  ARRAY, NO) 


*  MAXNUM  IS  THE  NUMBER  IN  THE  VECTOR  TO  BE 

*  SORTED 

INTEGER  MAXNUM.OFFSET.SW, NO, LIMIT 
DIMENSION  NO(MAXNUM) 

REAL  ARRAY(MAXNUM),TEMP,TEMP2 
LOGICAL  SWITCH 

* 

PRINT  *,  ’SUB  SHELL’  OFFSET=MAXNUM/2 
1  IF  (OFFSET  .GT.  0)  THEN 
LIMIT=MAXNUM-OFFSET 
SWITCH=  .FALSE. 


DO  15  J=l, LIMIT 

IF  (ARRAY(J)  .GE.  ARRAY(J+ OFFSET))  THEN 
TEMP  =  ARRAY(J) 

TEMP2  =  NO(J) 
ARRAY(J)=ARRAY(J+OFFSET) 
NO(J)=NO(J+OFFSET) 
ARRAY(J+OFFSET)=TEMP 
NO(J+OFFSET)=TEMP2 
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*  PERFORMS  A  SIMULATION  OF  EXPERIMENTS  USING  THE  * 

*  WEIBULL  MODEL 

* 

INTEGER  N 

REAL  ALPHA, BETA, MEAN,MAXSLK,RAND2(N), 
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+SLACK(N),RNDM 
MAXSLK  =  0.0 
PRINT  *,  ’SUB  NORMAL’ 

*-  CALCULATE  THE  STRENGTH  OF  THE  WEIBULL  FUNCTION, 
DO  10  I=1N 
RNDM  =  RAND2(I) 

S  L  A  C  K  (  I  )  = 

+EXP((LOG(-LOG(l-RNDM))+ALPHA*LOG(BETA))/ALPHA) 

IF  (SLACK(I)  .GT.  MAXSLK)  THEN 
MAXSLK  =  SLACK(I) 

END  IF 

10  CONTINUE 
END 

*  END 


SUBROUTINE 

+  UNIFM(RANGE1,RANGE2,N,RAND2,SLACK,MAXSLK) 


INTEGER  N 

REAL  RANGE,ME AN, RAND2(N),SLACK(N), MAXSLK 

* 

CONST=RANGE2-RANGE  1 

MAXSLK=0.0 

DO  10  I=1,N 

SLACK(I)=RAND2(I)*CONST  +  RANGE1 
IF  (MAXSLK  .LT.  SLACK(I))  THEN 
MAXSLK=SLACK(  I) 

END  IF 

10  CONTINUE 
END 


SUBROUTINE  WEIBUUALPHA, BETA, TOTAL, BUNDLE, P) 


*  PERFORMS  A  SIMULATION  OF  EXPERIMENTS  USING  THE 

*  WEIBULL  MODEL 

*  INITIALIZE 
REAL  ALPHA.BETA 
INTEGER  TOTAL.MAX 

REAL  P(TOTAL), BUND LE(TOTAL), RNDM 
PRINT  *  ’SUB  WEIBUL’ 

CALCULATE  THE  STRENGTH  OF  THE  WEIBULL  FUNCTION, 
DO  10  1=1,  TOTAL 
RNDM  =  BUNDLE(I) 

P(I)  = 

+  EXP((LOG{-LOG<l-RNDM))+ALPHA*LOG(BETA))/ALPHA) 

10  CONTINUE 
END 
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